From 06d6ecfd7f5b873f904cf58d1c862b38783454ec Mon Sep 17 00:00:00 2001 From: Jim Carr Date: Sun, 13 Dec 2020 13:57:56 -0500 Subject: [PATCH] Calculate -> Local moonrise and moonset. --- PALib.Tests/PAMoon.cs | 6 + PALib/PAMacros.cs | 719 ++++++++++++++++++++++++++++++++++++++++++ PALib/PAMoon.cs | 45 +++ README.md | 2 +- 4 files changed, 771 insertions(+), 1 deletion(-) diff --git a/PALib.Tests/PAMoon.cs b/PALib.Tests/PAMoon.cs index 5064772..8cb9800 100644 --- a/PALib.Tests/PAMoon.cs +++ b/PALib.Tests/PAMoon.cs @@ -40,5 +40,11 @@ public void MoonDistAngDiamHorParallax() { Assert.Equal((367964, 0, 32, 0, 59, 35.49), _paMoon.MoonDistAngDiamHorParallax(0, 0, 0, false, 0, 1, 9, 2003)); } + + [Fact] + public void MoonriseAndMoonset() + { + Assert.Equal((4, 21, 6, 3, 1986, 127.34, 13, 8, 6, 3, 1986, 234.05), _paMoon.MoonriseAndMoonset(6, 3, 1986, false, -5, -71.05, 42.3667)); + } } } \ No newline at end of file diff --git a/PALib/PAMacros.cs b/PALib/PAMacros.cs index 69422e0..f92b43d 100644 --- a/PALib/PAMacros.cs +++ b/PALib/PAMacros.cs @@ -4289,5 +4289,724 @@ public static double Sign(double numberToCheck) return signValue; } + + /// + /// Original macro name: UTDayAdjust + /// + public static double UTDayAdjust(double ut, double g1) + { + var returnValue = ut; + + if ((ut - g1) < -6.0) + returnValue = ut + 24.0; + + if ((ut - g1) > 6.0) + returnValue = ut - 24.0; + + return returnValue; + } + + /// + /// Original macro name: Fpart + /// + public static double FPart(double w) + { + return w - Lint(w); + } + + /// + /// Local time of moonrise. + /// + /// + /// Original macro name: MoonRiseLCT + /// + /// + /// + /// + /// + /// + /// + /// + /// hours + public static double MoonRiseLCT(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var lct6700result1 = MoonRiseLCT_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = lct6700result1.lu; + lct = lct6700result1.lct; + + if (lct == -99.0) + return lct; + + var la = lu; + + double x; + double ut; + var g1 = 0.0; + var gu = 0.0; + + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var lct6680result = MoonRiseLCT_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = lct6680result.lct; + dy1 = lct6680result.dy1; + mn1 = lct6680result.mn1; + yr1 = lct6680result.yr1; + gdy = lct6680result.gdy; + gmn = lct6680result.gmn; + gyr = lct6680result.gyr; + + var lct6700result2 = MoonRiseLCT_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + lu = lct6700result2.lu; + lct = lct6700result2.lct; + + if (lct == -99.0) + return lct; + + la = lu; + } + + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + + return lct; + } + + /// + /// Helper function for MoonRiseLCT + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonRiseLCT_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for MoonRiseLCT + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct) MoonRiseLCT_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeRise(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + if (!ERS(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat).Equals("OK")) + lct = -99.0; + + return (mm, bm, pm, dp, th, di, p, q, lu, lct); + } + + /// + /// Local date of moonrise. + /// + /// + /// Original macro names: MoonRiseLcDay, MoonRiseLcMonth, MoonRiseLcYear + /// + /// + /// Local date (day) + /// Local date (month) + /// Local date (year) + /// + public static (double dy1, int mn1, int yr1) MoonRiseLcDMY(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var lct6700result1 = MoonRiseLcDMY_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = lct6700result1.lu; + lct = lct6700result1.lct; + + if (lct == -99.0) + return (lct, (int)lct, (int)lct); + + var la = lu; + + double x; + double ut; + var g1 = 0.0; + var gu = 0.0; + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var lct6680result1 = MoonRiseLcDMY_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = lct6680result1.lct; + dy1 = lct6680result1.dy1; + mn1 = lct6680result1.mn1; + yr1 = lct6680result1.yr1; + gdy = lct6680result1.gdy; + gmn = lct6680result1.gmn; + gyr = lct6680result1.gyr; + + var lct6700result2 = MoonRiseLcDMY_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + + lu = lct6700result2.lu; + lct = lct6700result2.lct; + + if (lct == -99.0) + return (lct, (int)lct, (int)lct); + + la = lu; + } + + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + + return (dy1, mn1, yr1); + } + + /// + /// Helper function for MoonRiseLcDMY + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonRiseLcDMY_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for MoonRiseLcDMY + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct) MoonRiseLcDMY_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeRise(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + return (mm, bm, pm, dp, th, di, p, q, lu, lct); + } + + /// + /// Local azimuth of moonrise. + /// + /// + /// Original macro name: MoonRiseAz + /// + /// degrees + public static double MoonRiseAz(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var az6700result1 = MoonRiseAz_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = az6700result1.lu; + lct = az6700result1.lct; + double au; + + if (lct == -99.0) + return lct; + + var la = lu; + + double x; + double ut; + double g1; + var gu = 0.0; + var aa = 0.0; + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var az6680result1 = MoonRiseAz_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = az6680result1.lct; + dy1 = az6680result1.dy1; + mn1 = az6680result1.mn1; + yr1 = az6680result1.yr1; + gdy = az6680result1.gdy; + gmn = az6680result1.gmn; + gyr = az6680result1.gyr; + + var az6700result2 = MoonRiseAz_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + lu = az6700result2.lu; + lct = az6700result2.lct; + au = az6700result2.au; + + if (lct == -99.0) + return lct; + + la = lu; + aa = au; + } + + au = aa; + + return au; + } + + /// + /// Helper function for MoonRiseAz + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonRiseAz_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for MoonRiseAz + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct, double au) MoonRiseAz_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeRise(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + var au = RiseSetAzimuthRise(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + return (mm, bm, pm, dp, th, di, p, q, lu, lct, au); + } + + /// + /// Local time of moonset. + /// + /// + /// Original macro name: MoonSetLCT + /// + /// hours + public static double MoonSetLCT(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var lct6700result1 = MoonSetLCT_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = lct6700result1.lu; + lct = lct6700result1.lct; + + if (lct == -99.0) + return lct; + + var la = lu; + + double x; + double ut; + var g1 = 0.0; + var gu = 0.0; + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var lct6680result1 = MoonSetLCT_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = lct6680result1.lct; + dy1 = lct6680result1.dy1; + mn1 = lct6680result1.mn1; + yr1 = lct6680result1.yr1; + gdy = lct6680result1.gdy; + gmn = lct6680result1.gmn; + gyr = lct6680result1.gyr; + + var lct6700result2 = MoonSetLCT_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + lu = lct6700result2.lu; + lct = lct6700result2.lct; + + if (lct == -99.0) + return lct; + + la = lu; + } + + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + + return lct; + } + + /// + /// Helper function for MoonSetLCT + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonSetLCT_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for MoonSetLCT + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct) MoonSetLCT_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeSet(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + if (!ERS(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat).Equals("OK")) + lct = -99.0; + + return (mm, bm, pm, dp, th, di, p, q, lu, lct); + } + + /// + /// Local date of moonset. + /// + /// + /// Original macro names: MoonSetLcDay, MoonSetLcMonth, MoonSetLcYear + /// + /// + /// Local date (day) + /// Local date (month) + /// Local date (year) + /// + public static (double dy1, int mn1, int yr1) MoonSetLcDMY(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var dmy6700result1 = MoonSetLcDMY_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = dmy6700result1.lu; + lct = dmy6700result1.lct; + + if (lct == -99.0) + return (lct, (int)lct, (int)lct); + + var la = lu; + + double x; + double ut; + var g1 = 0.0; + var gu = 0.0; + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var dmy6680result1 = MoonSetLcDMY_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = dmy6680result1.lct; + dy1 = dmy6680result1.dy1; + mn1 = dmy6680result1.mn1; + yr1 = dmy6680result1.yr1; + gdy = dmy6680result1.gdy; + gmn = dmy6680result1.gmn; + gyr = dmy6680result1.gyr; + + var dmy6700result2 = MoonSetLcDMY_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + lu = dmy6700result2.lu; + lct = dmy6700result2.lct; + + if (lct == -99.0) + return (lct, (int)lct, (int)lct); + + la = lu; + } + + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + + return (dy1, mn1, yr1); + } + + /// + /// Helper function for MoonSetLcDMY + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonSetLcDMY_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for MoonSetLcDMY + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct) MoonSetLcDMY_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeSet(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + return (mm, bm, pm, dp, th, di, p, q, lu, lct); + } + + /// + /// Local azimuth of moonset. + /// + /// + /// Original macro name: MoonSetAz + /// + /// degrees + public static double MoonSetAz(double dy, int mn, int yr, int ds, int zc, double gLong, double gLat) + { + var gdy = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gmn = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var gyr = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var lct = 12.0; + var dy1 = dy; + var mn1 = mn; + var yr1 = yr; + + var az6700result1 = MoonSetAz_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + var lu = az6700result1.lu; + lct = az6700result1.lct; + + double au; + + if (lct == -99.0) + return lct; + + var la = lu; + + double x; + double ut; + double g1; + var gu = 0.0; + var aa = 0.0; + for (int k = 1; k < 9; k++) + { + x = LocalSiderealTimeToGreenwichSiderealTime(la, 0.0, 0.0, gLong); + ut = GreenwichSiderealTimeToUniversalTime(x, 0.0, 0.0, gdy, gmn, gyr); + + g1 = (k == 1) ? ut : gu; + + gu = ut; + ut = gu; + + var az6680result1 = MoonSetAz_L6680(x, ds, zc, gdy, gmn, gyr, g1, ut); + lct = az6680result1.lct; + dy1 = az6680result1.dy1; + mn1 = az6680result1.mn1; + yr1 = az6680result1.yr1; + gdy = az6680result1.gdy; + gmn = az6680result1.gmn; + gyr = az6680result1.gyr; + + var az6700result2 = MoonSetAz_L6700(lct, ds, zc, dy1, mn1, yr1, gdy, gmn, gyr, gLat); + lu = az6700result2.lu; + lct = az6700result2.lct; + au = az6700result2.au; + + if (lct == -99.0) + return lct; + + la = lu; + aa = au; + } + + au = aa; + + return au; + } + + /// + /// Helper function for moon_set_az + /// + public static (double ut, double lct, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr) MoonSetAz_L6680(double x, int ds, int zc, double gdy, int gmn, int gyr, double g1, double ut) + { + if (!EGstUt(x, 0.0, 0.0, gdy, gmn, gyr).Equals("OK")) + if (Math.Abs(g1 - ut) > 0.5) + ut = ut + 23.93447; + + ut = UTDayAdjust(ut, g1); + var lct = UniversalTimeToLocalCivilTime(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var dy1 = UniversalTime_LocalCivilDay(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var mn1 = UniversalTime_LocalCivilMonth(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + var yr1 = UniversalTime_LocalCivilYear(ut, 0.0, 0.0, ds, zc, gdy, gmn, gyr); + gdy = LocalCivilTimeGreenwichDay(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gmn = LocalCivilTimeGreenwichMonth(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + gyr = LocalCivilTimeGreenwichYear(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + ut = ut - 24.0 * (ut / 24.0).Floor(); + + return (ut, lct, dy1, mn1, yr1, gdy, gmn, gyr); + } + + /// + /// Helper function for moon_set_az + /// + public static (double mm, double bm, double pm, double dp, double th, double di, double p, double q, double lu, double lct, double au) MoonSetAz_L6700(double lct, int ds, int zc, double dy1, int mn1, int yr1, double gdy, int gmn, int gyr, double gLat) + { + var mm = MoonLong(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var bm = MoonLat(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1); + var pm = (MoonHP(lct, 0.0, 0.0, ds, zc, dy1, mn1, yr1)).ToRadians(); + var dp = NutatLong(gdy, gmn, gyr); + var th = 0.27249 * pm.Sine(); + var di = th + 0.0098902 - pm; + var p = DecimalDegreesToDegreeHours(EcRA(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr)); + var q = EcDec(mm + dp, 0.0, 0.0, bm, 0.0, 0.0, gdy, gmn, gyr); + var lu = RiseSetLocalSiderealTimeSet(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + var au = RiseSetAzimuthSet(p, 0.0, 0.0, q, 0.0, 0.0, Degrees(di), gLat); + + return (mm, bm, pm, dp, th, di, p, q, lu, lct, au); + } + } } \ No newline at end of file diff --git a/PALib/PAMoon.cs b/PALib/PAMoon.cs index da10216..42915d8 100644 --- a/PALib/PAMoon.cs +++ b/PALib/PAMoon.cs @@ -273,5 +273,50 @@ public class PAMoon return (earthMoonDist, angDiameterDeg, angDiameterMin, horParallaxDeg, horParallaxMin, horParallaxSec); } + + /// + /// Calculate date/time of local moonrise and moonset. + /// + /// + /// mrLTHour -- Moonrise, local time (hour part) + /// mrLTMin -- Moonrise, local time (minutes part) + /// mrLocalDateDay -- Moonrise, local date (day) + /// mrLocalDateMonth -- Moonrise, local date (month) + /// mrLocalDateYear -- Moonrise, local date (year) + /// mrAzimuthDeg -- Moonrise, azimuth (degrees) + /// msLTHour -- Moonset, local time (hour part) + /// msLTMin -- Moonset, local time (minutes part) + /// msLocalDateDay -- Moonset, local date (day) + /// msLocalDateMonth -- Moonset, local date (month) + /// msLocalDateYear -- Moonset, local date (year) + /// msAzimuthDeg -- Moonset, azimuth (degrees) + /// + public (double mrLTHour, double mrLTMin, double mrLocalDateDay, int mrLocalDateMonth, int mrLocalDateYear, double mrAzimuthDeg, double msLTHour, double msLTMin, double msLocalDateDay, int msLocalDateMonth, int msLocalDateYear, double msAzimuthDeg) MoonriseAndMoonset(double localDateDay, int localDateMonth, int localDateYear, bool isDaylightSaving, int zoneCorrectionHours, double geogLongDeg, double geogLatDeg) + { + var daylightSaving = (isDaylightSaving) ? 1 : 0; + + var localTimeOfMoonriseHours = PAMacros.MoonRiseLCT(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + var moonRiseLCResult = PAMacros.MoonRiseLcDMY(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + var localAzimuthDeg1 = PAMacros.MoonRiseAz(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + + var localTimeOfMoonsetHours = PAMacros.MoonSetLCT(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + var moonSetLCResult = PAMacros.MoonSetLcDMY(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + var localAzimuthDeg2 = PAMacros.MoonSetAz(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongDeg, geogLatDeg); + + var mrLTHour = PAMacros.DecimalHoursHour(localTimeOfMoonriseHours + 0.008333); + var mrLTMin = PAMacros.DecimalHoursMinute(localTimeOfMoonriseHours + 0.008333); + var mrLocalDateDay = moonRiseLCResult.dy1; + var mrLocalDateMonth = moonRiseLCResult.mn1; + var mrLocalDateYear = moonRiseLCResult.yr1; + var mrAzimuthDeg = Math.Round(localAzimuthDeg1, 2); + var msLTHour = PAMacros.DecimalHoursHour(localTimeOfMoonsetHours + 0.008333); + var msLTMin = PAMacros.DecimalHoursMinute(localTimeOfMoonsetHours + 0.008333); + var msLocalDateDay = moonSetLCResult.dy1; + var msLocalDateMonth = moonSetLCResult.mn1; + var msLocalDateYear = moonSetLCResult.yr1; + var msAzimuthDeg = Math.Round(localAzimuthDeg2, 2); + + return (mrLTHour, mrLTMin, mrLocalDateDay, mrLocalDateMonth, mrLocalDateYear, mrAzimuthDeg, msLTHour, msLTMin, msLocalDateDay, msLocalDateMonth, msLocalDateYear, msAzimuthDeg); + } } } \ No newline at end of file diff --git a/README.md b/README.md index a6ed944..9997880 100644 --- a/README.md +++ b/README.md @@ -60,7 +60,7 @@ If you're interested in this topic, please buy the book! It provides far more de - [x] Calculate -> Moon phase and position angle of bright limb - [x] Calculate -> Times of new Moon and full Moon - [x] Calculate -> Moon's distance, angular diameter, and horizontal parallax -- [ ] Calculate -> Local moonrise and moonset +- [x] Calculate -> Local moonrise and moonset ### Eclipses