From b9a4824fe387552ba2481c42caca2cb344d83b30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Richard=20K=C3=B6rber?= Date: Sun, 3 Mar 2024 09:48:27 +0100 Subject: [PATCH] Fix strange results on topocentric mode --- .../commons/suncalc/MoonIllumination.java | 28 ++++++++++--------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/main/java/org/shredzone/commons/suncalc/MoonIllumination.java b/src/main/java/org/shredzone/commons/suncalc/MoonIllumination.java index b9dbf39..f8a322b 100644 --- a/src/main/java/org/shredzone/commons/suncalc/MoonIllumination.java +++ b/src/main/java/org/shredzone/commons/suncalc/MoonIllumination.java @@ -93,14 +93,8 @@ public Parameters geocentric() { @Override public MoonIllumination execute() { JulianDate t = getJulianDate(); - Vector s, m; - if (hasLocation()) { - s = Sun.positionTopocentric(t, getLatitudeRad(), getLongitudeRad(), getElevation()); - m = Moon.positionTopocentric(t, getLatitudeRad(), getLongitudeRad(), getElevation()); - } else { - s = Sun.position(t); - m = Moon.position(t); - } + Vector s = Sun.position(t); + Vector m = Moon.position(t); double phi = PI - acos(m.dot(s) / (m.getR() * s.getR())); Vector sunMoon = m.cross(s); @@ -109,12 +103,20 @@ public MoonIllumination execute() { sin(s.getTheta()) * cos(m.getTheta()) - cos(s.getTheta()) * sin(m.getTheta()) * cos(s.getPhi() - m.getPhi()) ); - double r = m.subtract(s).norm(); - double re = s.norm(); - double d = m.norm(); - double elongation = acos((d*d + re*re - r*r) / (2.0*d*re)); + Vector sTopo, mTopo; + if (hasLocation()) { + sTopo = Sun.positionTopocentric(t, getLatitudeRad(), getLongitudeRad(), getElevation()); + mTopo = Moon.positionTopocentric(t, getLatitudeRad(), getLongitudeRad(), getElevation()); + } else { + sTopo = s; + mTopo = m; + } - double moonRadius = Moon.angularRadius(m.getR()); + double r = mTopo.subtract(sTopo).norm(); + double re = sTopo.norm(); + double d = mTopo.norm(); + double elongation = acos((d*d + re*re - r*r) / (2.0*d*re)); + double moonRadius = Moon.angularRadius(mTopo.getR()); double crescentWidth = moonRadius * (1 - cos(elongation)); return new MoonIllumination(