diff --git a/src/plugins/layers/useGrayLine.js b/src/plugins/layers/useGrayLine.js index cf7b17e..c83bc4d 100644 --- a/src/plugins/layers/useGrayLine.js +++ b/src/plugins/layers/useGrayLine.js @@ -103,42 +103,70 @@ function calculateSolarAltitude(date, latitude, longitude) { function generateTerminatorLine(date, solarAltitude = 0, numPoints = 360) { const points = []; const { declination } = calculateSolarPosition(date); + const decRad = declination * Math.PI / 180; + const altRad = solarAltitude * Math.PI / 180; + // For each longitude, calculate the latitude where the sun is at the specified altitude for (let i = 0; i <= numPoints; i++) { const lon = (i / numPoints) * 360 - 180; - - // Calculate latitude where sun is at specified altitude const hourAngle = calculateHourAngle(date, lon); const haRad = hourAngle * Math.PI / 180; - const decRad = declination * Math.PI / 180; - const altRad = solarAltitude * Math.PI / 180; - // Solve for latitude using solar altitude equation + // Use the solar altitude equation to solve for latitude + // sin(altitude) = sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(HA) + // Rearranging: sin(altitude) - sin(lat) * sin(dec) = cos(lat) * cos(dec) * cos(HA) + + // For terminator (altitude = 0), the equation simplifies + // We need to solve: tan(lat) = -tan(dec) / cos(HA) + const cosHA = Math.cos(haRad); const sinDec = Math.sin(decRad); const cosDec = Math.cos(decRad); const sinAlt = Math.sin(altRad); - // lat = arcsin((sin(alt) - sin(dec) * sin(lat)) / (cos(dec) * cos(lat))) - // Simplified: lat where sun altitude equals target - const numerator = sinAlt - sinDec * Math.sin(0); // approximation - const denominator = cosDec * cosHA; + // Solve using the quadratic formula or direct calculation + // sin(lat) = (sin(alt) - cos(lat) * cos(dec) * cos(HA)) / sin(dec) + + // Better approach: use atan2 for proper terminator calculation + // The terminator latitude for a given longitude is: + // lat = atan(-cos(HA) / tan(dec)) when dec != 0 let lat; - if (Math.abs(denominator) < 0.001) { - lat = declination > 0 ? 90 : -90; + + if (Math.abs(declination) < 0.01) { + // Near equinox: terminator is nearly straight along equator + lat = 0; } else { - const tanLat = (sinAlt - sinDec * Math.sin(declination * Math.PI / 180)) / denominator; - lat = Math.atan(tanLat) * 180 / Math.PI; + // Standard case: calculate terminator latitude + // Formula: cos(lat) * cos(dec) * cos(HA) = -sin(lat) * sin(dec) (for altitude = 0) + // This gives: tan(lat) = -cos(HA) / tan(dec) + + const tanDec = Math.tan(decRad); + if (Math.abs(tanDec) < 0.0001) { + lat = 0; + } else { + lat = Math.atan(-cosHA / tanDec) * 180 / Math.PI; + } - // More accurate calculation - const cosLat = Math.cos(lat * Math.PI / 180); - const sinLat = Math.sin(lat * Math.PI / 180); - const recalc = sinLat * sinDec + cosLat * cosDec * cosHA; - lat = Math.asin(Math.max(-1, Math.min(1, recalc))) * 180 / Math.PI; + // For twilight (altitude < 0), we need to adjust + if (solarAltitude !== 0) { + // Use iterative solution for twilight calculations + // cos(lat) * cos(dec) * cos(HA) + sin(lat) * sin(dec) = sin(alt) + + // Newton-Raphson iteration to solve for latitude + let testLat = lat * Math.PI / 180; + for (let iter = 0; iter < 5; iter++) { + const f = Math.sin(testLat) * sinDec + Math.cos(testLat) * cosDec * cosHA - sinAlt; + const fPrime = Math.cos(testLat) * sinDec - Math.sin(testLat) * cosDec * cosHA; + if (Math.abs(fPrime) > 0.0001) { + testLat = testLat - f / fPrime; + } + } + lat = testLat * 180 / Math.PI; + } } - // Clamp latitude + // Clamp latitude to valid range lat = Math.max(-90, Math.min(90, lat)); if (isFinite(lat) && isFinite(lon)) {