Commit: 9fb221ae7e342c149d465d994810a94fcc6730c0 Author: Vi Grey Date: 2023-10-11 10:23 UTC Summary: Better eclipse tracking and fractional seconds - Fix "true" and "false" of eclipse outputs to not have quotes - Calculate inner and outer shared tangents of Sun and Earth for eclipse detection - Use fractional seconds for better accuracy README.md | 58 +++++++++++++++++++++++++++++----------------------------- syzygy.hs | 61 +++++++++++++++++++++++++++++++++++-------------------------- 2 files changed, 64 insertions(+), 55 deletions(-) diff --git a/README.md b/README.md index 56116ef..f49dd4e 100644 --- a/README.md +++ b/README.md @@ -59,57 +59,57 @@ In cases where **\** and **\** are not specified, syzygy will use 0 ## Example Usage ```sh -syzygy 1681960707 +syzygy 42.9913021 -87.8838436 1681960707.141592 ``` The above command will provide the following output ``` { - "moonIllumination": 3.2849687072900036e-5, - "moonPhase": -0.4906395450400396, - "sunAngle": 50.80166346578888, - "moonAngle": 50.265302577732506, + "moonIllumination": 3.2847973790328666e-5, + "moonPhase": -0.4906185698952754, + "sunAngle": 64.38584608284015, + "moonAngle": 62.869180780600075, "sunGeoditicLocation": { - "latitude": 11.400793776368719, - "longitude": 130.14335777781847 + "latitude": 11.400794340426048, + "longitude": 130.14276760558175 }, "moonGeoditicLocation": { - "latitude": 10.822721967150077, - "longitude": 129.8314815597577 + "latitude": 10.822731885063867, + "longitude": 129.83091061511882 }, "earthGeoditicLocation": { - "latitude": 0.0, - "longitude": 0.0 + "latitude": 42.9913021, + "longitude": -87.8838436 }, "sunECEF": { - "x": -9.496465728561985e7, - "y": 1.1260110346287003e8, - "z": 2.9703008309426412e7 + "x": -9.496349729293042e7, + "y": 1.1260208146263802e8, + "z": 2.9703009772917412e7 }, "moonECEF": { - "x": -236404.3977138535, - "y": 283424.5644229226, - "z": 70548.5495178427 + "x": -236401.5695450799, + "y": 283426.9155068214, + "z": 70548.61457743934 }, "earthECEF": { - "x": 6378.137, - "y": 0.0, - "z": 0.0 + "x": 172.5385448589644, + "y": -4669.42563005196, + "z": 4326.7950220101175 }, "equationOfTime": { - "eccentricity": -1.8414647685022605, - "obliquity": 2.0806876889020387 + "eccentricity": -1.8414647539088742, + "obliquity": 2.0806877612331327 }, "eclipse": { - "local": "false", - "global": "true" + "local": false, + "global": true }, - "sunEarthPhase": -130.1433577778185, - "sunEclipticLongitude": 29.798841360580102, - "moonEclipticLongitude": 29.308201815540087, - "sunEarthMoonAngle": 0.6551387342734636, - "time": 1.681960707e9 + "sunEarthPhase": 141.97338879441824, + "sunEclipticLongitude": 29.79884296235347, + "moonEclipticLongitude": 29.308224392458214, + "sunEarthMoonAngle": 0.6551216493869961, + "time": 1.681960707141592e9 } ``` diff --git a/syzygy.hs b/syzygy.hs index 123beef..81d9540 100644 --- a/syzygy.hs +++ b/syzygy.hs @@ -875,13 +875,19 @@ earthLatLonArgs args jdeArg :: [String] -> POSIXTime -> (Double, Double) jdeArg args curTime - | length args > 2 = (jdeOfUnixTimestamp (fromIntegral (floor (read (args !! 2) :: Double))), fromIntegral (floor (read (args !! 2) :: Double))) - | length args == 1 = (jdeOfUnixTimestamp (fromIntegral (floor (read (args !! 0) :: Double))), fromIntegral (floor (read (args !! 0) :: Double))) - | otherwise = (jdeOfUnixTimestamp (fromIntegral (floor curTime)), fromIntegral (floor curTime)) + | length args > 2 = (jdeOfUnixTimestamp (read (args !! 2) :: Double), read (args !! 2) :: Double) + | length args == 1 = (jdeOfUnixTimestamp (read (args !! 0) :: Double), read (args !! 0) :: Double) + | otherwise = (jdeOfUnixTimestamp (fromIntegral (floor (curTime * 1000)) / 1000), fromIntegral (floor (curTime * 1000)) / 1000) --- Determines if point C is to the "left" of vector A->B -isLeftSide :: Double -> Double -> Double -> Double -> Double -> Double -> Bool -isLeftSide aX aY bX bY cX cY = (bX-aX)*(cY-aY)-(bY-aY)*(cX-aX) < 0 +-- Determine if (cX,cY) is left of the line defined by slope m and +-- y-intercept b +isLeftSide :: Double -> Double -> Double -> Double -> Bool +isLeftSide m b cX cY = (cX - cY / m) < ((-b) / m) + +-- Determine if the line defined by slope m and y-intercept b goes +-- through circle with radius of cR at (cX,cY) +isCircleCollision :: Double -> Double -> Double -> Double -> Double -> Bool +isCircleCollision m b cX cY cR = abs ((-m) * cX + cY - b) / sqrt (m**2 + 1) < cR -- Determine if c blocks a from b if c is in front of b or if b blocks -- a from c if c is behind b @@ -889,26 +895,30 @@ isLeftSide aX aY bX bY cX cY = (bX-aX)*(cY-aY)-(bY-aY)*(cX-aX) < 0 -- a, b, and c. isEclipse :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Bool isEclipse aX aY aZ aR bX bY bZ bR cX cY cZ cR - | cY0 > 0 && (isLeftSide aX2 aY2 bR 0 cX1 cY1) = True - | cY0 < 0 && (isLeftSide aX1 aY1 bR 0 cX2 cY2) = True + | cY0 > 0 && ((isCircleCollision outerM outerB cX0 cY0 cR) || (isLeftSide outerM outerB cX0 cY0)) = True + | cY0 < 0 && ((isCircleCollision innerM innerB cX0 cY0 cR) || (isLeftSide innerM innerB cX0 cY0)) = True | otherwise = False where lenAB = sqrt ((aX-bX)**2 + (aY-bY)**2 + (aZ-bZ)**2) - lenAC = sqrt ((aX-cX)**2 + (aY-cY)**2 + (aZ-cZ)**2) - lenBC = sqrt ((cX-bX)**2 + (cY-bY)**2 + (cZ-bZ)**2) - angleB = pi/2 - acos ((lenAB**2 + lenBC**2 - lenAC**2) / (2 * lenAB * lenBC)) - cX0 = cos (angleB) * lenBC - cY0 = sin (angleB) * lenBC - angleNewA = pi + (atan2 lenAB (-bR)) - angleNewC = pi + (atan2 cY0 (cX0 - bR)) - aX1 = cos (angleNewA - pi/2) * aR - aY1 = lenAB + sin (angleNewA - pi/2) * aR - aX2 = cos (angleNewA + pi/2) * aR - aY2 = lenAB + sin (angleNewA + pi/2) * aR - cX1 = cX0 + cos (angleNewC - pi/2) * cR - cY1 = cY0 + sin (angleNewC - pi/2) * cR - cX2 = cX0 + cos (angleNewC + pi/2) * cR - cY2 = cY0 + sin (angleNewC + pi/2) * cR + lenBC = sqrt ((bX-cX)**2 + (bY-cY)**2 + (bZ-cZ)**2) + lenCA = sqrt ((cX-aX)**2 + (cY-aY)**2 + (cZ-aZ)**2) + angleABC = pi/2 - acos ((lenAB**2 + lenBC**2 - lenCA**2) / (2 * lenAB * lenBC)) + cX0 = cos (angleABC) * lenBC + cY0 = sin (angleABC) * lenBC + outerTheta = pi * 2 - asin ((aR - bR) / lenAB) + innerTheta = (atan2 (-lenAB) 0) + asin ((aR + bR) / lenAB) - pi / 2 + aXO = cos (outerTheta) * aR + aYO = lenAB + sin (outerTheta) * aR + bXO = cos (outerTheta) * bR + bYO = sin (outerTheta) * bR + aXI = cos (innerTheta) * aR + aYI = lenAB + sin (innerTheta) * aR + bXI = cos (innerTheta) * bR + bYI = sin (innerTheta) * bR + outerM = (bYO - aYO) / (bXO - aXO) + outerB = aYO - outerM * aXO + innerM = (bYI - aYI) / (bXI - aXI) + innerB = aYI - innerM * aXI main = do setLocaleEncoding utf8 @@ -939,7 +949,6 @@ main = do let sunLamdaPrime = sunLonInit - t * 1.397 - t ** 2 * 0.00031 let dO = arcDegToDecDeg 0 0 (-0.09033) let dB = arcDegToDecDeg 0 0 0.03916 * (cos (degreesToRadians sunLamdaPrime) - sin (degreesToRadians sunLamdaPrime)) - let sunApparentLon = mod' (sunLonInit + dY + dO + (-0.005775518) * sunRadiusVector * dL) 360 let sunApparentLat = mod' (sunLatInit + dB) 360 let sunDist = auToKM sunRadiusVector @@ -1007,8 +1016,8 @@ main = do putStrLn (" \"obliquity\": " ++ show eotO) putStrLn (" },") putStrLn (" \"eclipse\": {") - putStrLn (" \"local\": " ++ show localEclipse ++ ",") - putStrLn (" \"global\": " ++ show globalEclipse) + putStrLn (" \"local\": " ++ localEclipse ++ ",") + putStrLn (" \"global\": " ++ globalEclipse) putStrLn (" },") putStrLn (" \"sunEarthPhase\": " ++ show (radiansToDegrees (degreesToRadians (earthLon - radiansToDegrees sunLon))) ++ ",") putStrLn (" \"sunEclipticLongitude\": " ++ show sunApparentLon ++ ",")