diff --git a/src/plugins/layers/useGrayLine.js b/src/plugins/layers/useGrayLine.js index 51df3b5..6ba37c9 100644 --- a/src/plugins/layers/useGrayLine.js +++ b/src/plugins/layers/useGrayLine.js @@ -100,6 +100,39 @@ function calculateSolarAltitude(date, latitude, longitude) { return altitude; } +// Split polyline at date line to avoid lines cutting across the map +function splitAtDateLine(points) { + if (points.length < 2) return [points]; + + const segments = []; + let currentSegment = [points[0]]; + + for (let i = 1; i < points.length; i++) { + const prev = points[i - 1]; + const curr = points[i]; + const prevLon = prev[1]; + const currLon = curr[1]; + const lonDiff = Math.abs(currLon - prevLon); + + // If longitude jumps more than 180°, we've crossed the date line + if (lonDiff > 180) { + // Finish current segment + segments.push(currentSegment); + // Start new segment + currentSegment = [curr]; + } else { + currentSegment.push(curr); + } + } + + // Add final segment + if (currentSegment.length > 0) { + segments.push(currentSegment); + } + + return segments.filter(seg => seg.length >= 2); +} + // Generate terminator line for a specific solar altitude function generateTerminatorLine(date, solarAltitude = 0, numPoints = 360) { const points = []; @@ -113,64 +146,78 @@ function generateTerminatorLine(date, solarAltitude = 0, numPoints = 360) { const hourAngle = calculateHourAngle(date, lon); const haRad = hourAngle * Math.PI / 180; - // 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); - // 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; + // Check if solution exists (sun can reach this altitude at this longitude) + // For terminator and twilight, check if |cos(HA) * cos(dec)| <= 1 - sin(alt) * sin(dec) + const testValue = (sinAlt - sinDec * sinDec) / (cosDec * cosDec * cosHA * cosHA); + if (Math.abs(declination) < 0.01) { // Near equinox: terminator is nearly straight along equator lat = 0; + } else if (Math.abs(cosDec) < 0.001) { + // Near solstice: sun is directly over tropic, skip this point + continue; } else { // 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; - } - // 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) + if (solarAltitude === 0) { + // Terminator (sunrise/sunset line) + // Formula: tan(lat) = -cos(HA) / tan(dec) + if (Math.abs(tanDec) > 0.0001) { + lat = Math.atan(-cosHA / tanDec) * 180 / Math.PI; + } else { + lat = 0; + } + } else { + // Twilight zones (negative solar altitude) + // Use Newton-Raphson iteration to solve for latitude + // Equation: sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(HA) = sin(alt) - // Newton-Raphson iteration to solve for latitude - let testLat = lat * Math.PI / 180; - for (let iter = 0; iter < 5; iter++) { + // Initial guess based on terminator + let testLat = Math.atan(-cosHA / tanDec); + + // Iterate to find solution + let converged = false; + for (let iter = 0; iter < 10; 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(f) < 0.0001) { + converged = true; + break; + } + if (Math.abs(fPrime) > 0.0001) { testLat = testLat - f / fPrime; + } else { + break; } + + // Constrain to valid latitude range during iteration + testLat = Math.max(-Math.PI/2, Math.min(Math.PI/2, testLat)); + } + + // Only use the point if iteration converged + if (!converged) { + continue; } + lat = testLat * 180 / Math.PI; } } - // Clamp latitude to valid range - lat = Math.max(-90, Math.min(90, lat)); + // Strict clamping to valid latitude range + lat = Math.max(-85, Math.min(85, lat)); - if (isFinite(lat) && isFinite(lon)) { + // Only add point if it's valid and not at extreme latitude + if (isFinite(lat) && isFinite(lon) && Math.abs(lat) < 85) { points.push([lat, lon]); } } @@ -488,104 +535,123 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) { // Main terminator (solar altitude = 0°) const terminator = generateTerminatorLine(currentTime, 0, 360); - const terminatorLine = L.polyline(terminator, { - color: '#ff6600', - weight: 3, - opacity: opacity * 0.8, - dashArray: '10, 5' + const terminatorSegments = splitAtDateLine(terminator); + + terminatorSegments.forEach(segment => { + const terminatorLine = L.polyline(segment, { + color: '#ff6600', + weight: 3, + opacity: opacity * 0.8, + dashArray: '10, 5' + }); + terminatorLine.bindPopup(` +