fix(grayline): Fix date line wrapping and Arctic line artifacts

TWO ISSUES FIXED: 1. Lines stopping at ±180° (international date line): - Added splitAtDateLine() function to detect date line crossings - Split polylines into segments when longitude jumps > 180° - Each segment renders separately, avoiding lines cutting across map 2. Random lines to Arctic (twilight zones): - Improved Newton-Raphson convergence checking - Added convergence validation before using calculated points - Stricter latitude clamping to ±85° to avoid polar discontinuities - Skip points near solstices where sun is directly over tropics - Filter out points with extreme latitudes (>85°) that cause artifacts - Validate enhanced DX zone has enough points before creating polygon Technical improvements: - More robust iterative solver for twilight zone calculations - Better handling of edge cases (equinox, solstice, polar regions) - Increased Newton-Raphson iterations from 5 to 10 for better accuracy - Added constraints during iteration to keep latitude in valid range Result: Clean, continuous lines across all longitudes including date line, no erratic lines near poles or Arctic regions.
1 day ago · ffd8a8ce4c
parent 47eb3ef107
commit ffd8a8ce4c
1 changed files with 174 additions and 108 deletions
--- 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 {
+      
+      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)
        
-      // 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)
+        // Initial guess based on terminator
+        let testLat = Math.atan(-cosHA / tanDec);
        
-        // Newton-Raphson iteration to solve for latitude
-        let testLat = lat * Math.PI / 180;
-        for (let iter = 0; iter < 5; iter++) {
+        // 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,7 +535,10 @@ 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, {
+    const terminatorSegments = splitAtDateLine(terminator);
+    
+    terminatorSegments.forEach(segment => {
+      const terminatorLine = L.polyline(segment, {
        color: '#ff6600',
        weight: 3,
        opacity: opacity * 0.8,
@ -504,12 +554,15 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) {
      `);
      terminatorLine.addTo(map);
      newLayers.push(terminatorLine);
+    });
    
    // Enhanced DX zone (±5° from terminator)
    if (showEnhancedZone) {
      const enhancedUpper = generateTerminatorLine(currentTime, 5, 360);
      const enhancedLower = generateTerminatorLine(currentTime, -5, 360);
      
+      // Only create polygon if we have valid points
+      if (enhancedUpper.length > 2 && enhancedLower.length > 2) {
        // Create polygon for enhanced zone
        const enhancedZone = [...enhancedUpper, ...enhancedLower.reverse()];
        const enhancedPoly = L.polygon(enhancedZone, {
@ -530,12 +583,16 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) {
        enhancedPoly.addTo(map);
        newLayers.push(enhancedPoly);
      }
+    }
    
    // Twilight zones
    if (showTwilight) {
      // Civil twilight (sun altitude -6°)
      const civilTwilight = generateTerminatorLine(currentTime, -6, 360);
-      const civilLine = L.polyline(civilTwilight, {
+      const civilSegments = splitAtDateLine(civilTwilight);
+      
+      civilSegments.forEach(segment => {
+        const civilLine = L.polyline(segment, {
          color: '#4488ff',
          weight: 2,
          opacity: twilightOpacity * 0.6,
@ -550,10 +607,14 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) {
        `);
        civilLine.addTo(map);
        newLayers.push(civilLine);
+      });
      
      // Nautical twilight (sun altitude -12°)
      const nauticalTwilight = generateTerminatorLine(currentTime, -12, 360);
-      const nauticalLine = L.polyline(nauticalTwilight, {
+      const nauticalSegments = splitAtDateLine(nauticalTwilight);
+      
+      nauticalSegments.forEach(segment => {
+        const nauticalLine = L.polyline(segment, {
          color: '#6666ff',
          weight: 1.5,
          opacity: twilightOpacity * 0.4,
@ -568,10 +629,14 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) {
        `);
        nauticalLine.addTo(map);
        newLayers.push(nauticalLine);
+      });
      
      // Astronomical twilight (sun altitude -18°)
      const astroTwilight = generateTerminatorLine(currentTime, -18, 360);
-      const astroLine = L.polyline(astroTwilight, {
+      const astroSegments = splitAtDateLine(astroTwilight);
+      
+      astroSegments.forEach(segment => {
+        const astroLine = L.polyline(segment, {
          color: '#8888ff',
          weight: 1,
          opacity: twilightOpacity * 0.3,
@ -586,6 +651,7 @@ export function useLayer({ enabled = false, opacity = 0.5, map = null }) {
        `);
        astroLine.addTo(map);
        newLayers.push(astroLine);
+      });
    }
    
    setLayers(newLayers);