From 3b44b2e27730283be92215934b53a82232aef174 Mon Sep 17 00:00:00 2001 From: Jim Carr Date: Wed, 16 Dec 2020 19:17:39 -0500 Subject: [PATCH] Calculate -> Solar eclipse occurrence and circumstances. --- PALib.Tests/PAEclipses.cs | 12 + PALib/PAEclipses.cs | 86 ++++++ PALib/PAMacros.cs | 578 ++++++++++++++++++++++++++++++++++++++ README.md | 2 +- 4 files changed, 677 insertions(+), 1 deletion(-) diff --git a/PALib.Tests/PAEclipses.cs b/PALib.Tests/PAEclipses.cs index f20d8e3..4902bbb 100644 --- a/PALib.Tests/PAEclipses.cs +++ b/PALib.Tests/PAEclipses.cs @@ -22,5 +22,17 @@ public void LunarEclipseCircumstances() { Assert.Equal((4, 4, 2015, 9, 0, 10, 16, 11, 55, 12, 1, 12, 7, 13, 46, 15, 1, 1.01), _paEclipses.LunarEclipseCircumstances(1, 4, 2015, false, 10)); } + + [Fact] + public void SolarEclipseOccurrence() + { + Assert.Equal(("Solar eclipse certain", 20, 3, 2015), _paEclipses.SolarEclipseOccurrence(1, 4, 2015, false, 0)); + } + + [Fact] + public void SolarEclipseCircumstances() + { + Assert.Equal((20, 3, 2015, 8, 55, 9, 57, 10, 58, 1.016), _paEclipses.SolarEclipseCircumstances(20, 3, 2015, false, 0, 0, 68.65)); + } } } \ No newline at end of file diff --git a/PALib/PAEclipses.cs b/PALib/PAEclipses.cs index 038bb1b..e0472d7 100644 --- a/PALib/PAEclipses.cs +++ b/PALib/PAEclipses.cs @@ -122,5 +122,91 @@ public class PAEclipses return (lunarEclipseCertainDateDay, lunarEclipseCertainDateMonth, lunarEclipseCertainDateYear, utStartPenPhaseHour, utStartPenPhaseMinutes, utStartUmbralPhaseHour, utStartUmbralPhaseMinutes, utStartTotalPhaseHour, utStartTotalPhaseMinutes, utMidEclipseHour, utMidEclipseMinutes, utEndTotalPhaseHour, utEndTotalPhaseMinutes, utEndUmbralPhaseHour, utEndUmbralPhaseMinutes, utEndPenPhaseHour, utEndPenPhaseMinutes, eclipseMagnitude); } + + /// + /// Determine if a solar eclipse is likely to occur. + /// + /// + /// status -- One of "Solar eclipse certain", "Solar eclipse possible", or "No solar eclipse". + /// event_date_day -- Date of eclipse event (day). + /// event_date_month -- Date of eclipse event (month). + /// event_date_year -- Date of eclipse event (year). + /// + public (string status, double eventDateDay, int eventDateMonth, int eventDateYear) SolarEclipseOccurrence(double localDateDay, int localDateMonth, int localDateYear, bool isDaylightSaving, int zoneCorrectionHours) + { + var daylightSaving = (isDaylightSaving) ? 1 : 0; + + var julianDateOfNewMoon = PAMacros.NewMoon(daylightSaving, zoneCorrectionHours, localDateDay, localDateMonth, localDateYear); + var gDateOfNewMoonDay = PAMacros.JulianDateDay(julianDateOfNewMoon); + var integerDay = (gDateOfNewMoonDay).Floor(); + var gDateOfNewMoonMonth = PAMacros.JulianDateMonth(julianDateOfNewMoon); + var gDateOfNewMoonYear = PAMacros.JulianDateYear(julianDateOfNewMoon); + var utOfNewMoonHours = gDateOfNewMoonDay - integerDay; + + var localCivilDateDay = PAMacros.UniversalTime_LocalCivilDay(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + var localCivilDateMonth = PAMacros.UniversalTime_LocalCivilMonth(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + var localCivilDateYear = PAMacros.UniversalTime_LocalCivilYear(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + + var eclipseOccurrence = PAMacros.SolarEclipseOccurrence(daylightSaving, zoneCorrectionHours, localDateDay, localDateMonth, localDateYear); + + var status = eclipseOccurrence; + var eventDateDay = localCivilDateDay; + var eventDateMonth = localCivilDateMonth; + var eventDateYear = localCivilDateYear; + + return (status, eventDateDay, eventDateMonth, eventDateYear); + } + + /// + /// Calculate the circumstances of a solar eclipse. + /// + /// + /// solarEclipseCertainDateDay -- Solar eclipse date (day) + /// solarEclipseCertainDateMonth -- Solar eclipse date (month) + /// solarEclipseCertainDateYear -- Solar eclipse date (year) + /// utFirstContactHour -- First contact of shadow (hour) + /// utFirstContactMinutes -- First contact of shadow (minutes) + /// utMidEclipseHour -- Mid-eclipse (hour) + /// utMidEclipseMinutes -- Mid-eclipse (minutes) + /// utLastContactHour -- Last contact of shadow (hour) + /// utLastContactMinutes -- Last contact of shadow (minutes) + /// eclipseMagnitude -- Eclipse magnitude + /// + public (double solarEclipseCertainDateDay, int solarEclipseCertainDateMonth, int solarEclipseCertainDateYear, double utFirstContactHour, double utFirstContactMinutes, double utMidEclipseHour, double utMidEclipseMinutes, double utLastContactHour, double utLastContactMinutes, double eclipseMagnitude) SolarEclipseCircumstances(double localDateDay, int localDateMonth, int localDateYear, bool isDaylightSaving, int zoneCorrectionHours, double geogLongitudeDeg, double geogLatitudeDeg) + { + var daylightSaving = (isDaylightSaving) ? 1 : 0; + + var julianDateOfNewMoon = PAMacros.NewMoon(daylightSaving, zoneCorrectionHours, localDateDay, localDateMonth, localDateYear); + var gDateOfNewMoonDay = PAMacros.JulianDateDay(julianDateOfNewMoon); + var integerDay = (gDateOfNewMoonDay).Floor(); + var gDateOfNewMoonMonth = PAMacros.JulianDateMonth(julianDateOfNewMoon); + var gDateOfNewMoonYear = PAMacros.JulianDateYear(julianDateOfNewMoon); + var utOfNewMoonHours = gDateOfNewMoonDay - integerDay; + var localCivilDateDay = PAMacros.UniversalTime_LocalCivilDay(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + var localCivilDateMonth = PAMacros.UniversalTime_LocalCivilMonth(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + var localCivilDateYear = PAMacros.UniversalTime_LocalCivilYear(utOfNewMoonHours, 0.0, 0.0, daylightSaving, zoneCorrectionHours, integerDay, gDateOfNewMoonMonth, gDateOfNewMoonYear); + + var utMaxEclipse = PAMacros.UTMaxSolarEclipse(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongitudeDeg, geogLatitudeDeg); + var utFirstContact = PAMacros.UTFirstContactSolarEclipse(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongitudeDeg, geogLatitudeDeg); + var utLastContact = PAMacros.UTLastContactSolarEclipse(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongitudeDeg, geogLatitudeDeg); + var magnitude = PAMacros.MagSolarEclipse(localDateDay, localDateMonth, localDateYear, daylightSaving, zoneCorrectionHours, geogLongitudeDeg, geogLatitudeDeg); + + var solarEclipseCertainDateDay = localCivilDateDay; + var solarEclipseCertainDateMonth = localCivilDateMonth; + var solarEclipseCertainDateYear = localCivilDateYear; + + var utFirstContactHour = (utFirstContact == -99.0) ? -99.0 : PAMacros.DecimalHoursHour(utFirstContact + 0.008333); + var utFirstContactMinutes = (utFirstContact == -99.0) ? -99.0 : PAMacros.DecimalHoursMinute(utFirstContact + 0.008333); + + var utMidEclipseHour = (utMaxEclipse == -99.0) ? -99.0 : PAMacros.DecimalHoursHour(utMaxEclipse + 0.008333); + var utMidEclipseMinutes = (utMaxEclipse == -99.0) ? -99.0 : PAMacros.DecimalHoursMinute(utMaxEclipse + 0.008333); + + var utLastContactHour = (utLastContact == -99.0) ? -99.0 : PAMacros.DecimalHoursHour(utLastContact + 0.008333); + var utLastContactMinutes = (utLastContact == -99.0) ? -99.0 : PAMacros.DecimalHoursMinute(utLastContact + 0.008333); + + var eclipseMagnitude = (magnitude == -99.0) ? -99.0 : Math.Round(magnitude, 3); + + return (solarEclipseCertainDateDay, solarEclipseCertainDateMonth, solarEclipseCertainDateYear, utFirstContactHour, utFirstContactMinutes, utMidEclipseHour, utMidEclipseMinutes, utLastContactHour, utLastContactMinutes, eclipseMagnitude); + } } } \ No newline at end of file diff --git a/PALib/PAMacros.cs b/PALib/PAMacros.cs index 95cd6fc..45f9b70 100644 --- a/PALib/PAMacros.cs +++ b/PALib/PAMacros.cs @@ -4305,6 +4305,34 @@ public static double FPart(double w) return w - Lint(w); } + /// + /// Original macro name: EQElat + /// + public static double EQELat(double rah, double ram, double ras, double dd, double dm, double ds, double gd, int gm, int gy) + { + var a = (DegreeHoursToDecimalDegrees(HMStoDH(rah, ram, ras))).ToRadians(); + var b = (DegreesMinutesSecondsToDecimalDegrees(dd, dm, ds)).ToRadians(); + var c = (Obliq(gd, gm, gy)).ToRadians(); + var d = b.Sine() * c.Cosine() - b.Cosine() * c.Sine() * a.Sine(); + + return Degrees(d.ASine()); + } + + /// + /// Original macro name: EQElong + /// + public static double EQELong(double rah, double ram, double ras, double dd, double dm, double ds, double gd, int gm, int gy) + { + var a = (DegreeHoursToDecimalDegrees(HMStoDH(rah, ram, ras))).ToRadians(); + var b = (DegreesMinutesSecondsToDecimalDegrees(dd, dm, ds)).ToRadians(); + var c = (Obliq(gd, gm, gy)).ToRadians(); + var d = a.Sine() * c.Cosine() + b.Tangent() * c.Sine(); + var e = a.Cosine(); + var f = Degrees(d.AngleTangent2(e)); + + return f - 360.0 * (f / 360.0).Floor(); + } + /// /// Local time of moonrise. /// @@ -5793,5 +5821,555 @@ public static double MagLunarEclipse(double dy, int mn, int yr, int ds, int zc) return mg; } + + /// + /// Determine if a solar eclipse is likely to occur. + /// + /// + /// Original macro name: SEOccurrence + /// + /// + /// + /// + /// + /// + /// + public static string SolarEclipseOccurrence(int ds, int zc, double dy, int mn, int yr) + { + var d0 = LocalCivilTimeGreenwichDay(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var m0 = LocalCivilTimeGreenwichMonth(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + var y0 = LocalCivilTimeGreenwichYear(12.0, 0.0, 0.0, ds, zc, dy, mn, yr); + + var j0 = CivilDateToJulianDate(0.0, 1, y0); + var dj = CivilDateToJulianDate(d0, m0, y0); + var k = (y0 - 1900.0 + ((dj - j0) * 1.0 / 365.0)) * 12.3685; + k = Lint(k + 0.5); + var tn = k / 1236.85; + var tf = (k + 0.5) / 1236.85; + var t = tn; + var l6855result1 = SolarEclipseOccurrence_L6855(t, k); + var nb = l6855result1.f; + t = tf; + k = k + 0.5; + var l6855result2 = SolarEclipseOccurrence_L6855(t, k); + + var df = Math.Abs(nb - 3.141592654 * Lint(nb / 3.141592654)); + + if (df > 0.37) + df = 3.141592654 - df; + + var s = "Solar eclipse certain"; + if (df >= 0.242600766) + { + s = "Solar eclipse possible"; + if (df > 0.37) + s = "No solar eclipse"; + } + + return s; + } + + /// + /// Helper function for SolarEclipseOccurrence + /// + public static (double f, double dd, double e1, double b1, double a, double b) SolarEclipseOccurrence_L6855(double t, double k) + { + var t2 = t * t; + var e = 29.53 * k; + var c = 166.56 + (132.87 - 0.009173 * t) * t; + c = c.ToRadians(); + var b = 0.00058868 * k + (0.0001178 - 0.000000155 * t) * t2; + b = b + 0.00033 * c.Sine() + 0.75933; + var a = k / 12.36886; + var a1 = 359.2242 + 360.0 * FPart(a) - (0.0000333 + 0.00000347 * t) * t2; + var a2 = 306.0253 + 360.0 * FPart(k / 0.9330851); + a2 = a2 + (0.0107306 + 0.00001236 * t) * t2; + a = k / 0.9214926; + var f = 21.2964 + 360.0 * FPart(a) - (0.0016528 + 0.00000239 * t) * t2; + a1 = UnwindDeg(a1); + a2 = UnwindDeg(a2); + f = UnwindDeg(f); + a1 = a1.ToRadians(); + a2 = a2.ToRadians(); + f = f.ToRadians(); + + var dd = (0.1734 - 0.000393 * t) * (a1).Sine() + 0.0021 * (2.0 * a1).Sine(); + dd = dd - 0.4068 * (a2).Sine() + 0.0161 * (2.0 * a2).Sine() - 0.0004 * (3.0 * a2).Sine(); + dd = dd + 0.0104 * (2.0 * f).Sine() - 0.0051 * (a1 + a2).Sine(); + dd = dd - 0.0074 * (a1 - a2).Sine() + 0.0004 * (2.0 * f + a1).Sine(); + dd = dd - 0.0004 * (2.0 * f - a1).Sine() - 0.0006 * (2.0 * f + a2).Sine() + 0.001 * (2.0 * f - a2).Sine(); + dd = dd + 0.0005 * (a1 + 2.0 * a2).Sine(); + var e1 = e.Floor(); + b = b + dd + (e - e1); + var b1 = b.Floor(); + a = e1 + b1; + b = b - b1; + + return (f, dd, e1, b1, a, b); + } + + /// + /// Calculate time of maximum shadow for solar eclipse (UT) + /// + /// + /// Original macro name: UTMaxSolarEclipse + /// + public static double UTMaxSolarEclipse(double dy, int mn, int yr, int ds, int zc, double glong, double glat) + { + var tp = 2.0 * Math.PI; + + if (SolarEclipseOccurrence(ds, zc, dy, mn, yr).Equals("No solar eclipse")) + return -99.0; + + var dj = NewMoon(ds, zc, dy, mn, yr); + var gday = JulianDateDay(dj); + var gmonth = JulianDateMonth(dj); + var gyear = JulianDateYear(dj); + var igday = gday.Floor(); + var xi = gday - igday; + var utnm = xi * 24.0; + var ut = utnm - 1.0; + var ly = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var my = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var by = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hy = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + ut = utnm + 1.0; + var sb = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians() - ly; + var mz = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var bz = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hz = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + + if (sb < 0.0) + sb = sb + tp; + + var xh = utnm; + var x = my; + var y = by; + var tm = xh - 1.0; + var hp = hy; + var l7390result1 = UTMaxSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + my = l7390result1.p; + by = l7390result1.q; + x = mz; + y = bz; + tm = xh + 1.0; + hp = hz; + var l7390result2 = UTMaxSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + mz = l7390result2.p; + bz = l7390result2.q; + + var x0 = xh + 1.0 - (2.0 * bz / (bz - by)); + var dm = mz - my; + + if (dm < 0.0) + dm = dm + tp; + + var lj = (dm - sb) / 2.0; + var mr = my + (dm * (x0 - xh + 1.0) / 2.0); + ut = x0 - 0.13851852; + var rr = SunDist(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear); + var sr = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + sr = sr + (NutatLong(igday, gmonth, gyear) - 0.00569).ToRadians(); + x = sr; + y = 0.0; + tm = ut; + hp = 0.00004263452 / rr; + var l7390result3 = UTMaxSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + // let(_paa, _qaa, _xaa, _pbb, _qbb, _xbb, p, q) = + sr = l7390result3.p; + by = by - l7390result3.q; + bz = bz - l7390result3.q; + var p3 = 0.00004263; + var zh = (sr - mr) / lj; + var tc = x0 + zh; + var sh = (((bz - by) * (tc - xh - 1.0) / 2.0) + bz) / lj; + var s2 = sh * sh; + var z2 = zh * zh; + var ps = p3 / (rr * lj); + var z1 = (zh * z2 / (z2 + s2)) + x0; + var h0 = (hy + hz) / (2.0 * lj); + var rm = 0.272446 * h0; + var rn = 0.00465242 / (lj * rr); + var hd = h0 * 0.99834; + var _ru = (hd - rn + ps) * 1.02; + var _rp = (hd + rn + ps) * 1.02; + var pj = Math.Abs(sh * zh / (s2 + z2).SquareRoot()); + var r = rm + rn; + var dd = z1 - x0; + dd = dd * dd - ((z2 - (r * r)) * dd / zh); + + if (dd < 0.0) + return -99.0; + + var zd = dd.SquareRoot(); + + return z1; + } + + /// + /// Helper function for ut_max_solar_eclipse + /// + public static (double paa, double qaa, double xaa, double pbb, double qbb, double xbb, double p, double q) UTMaxSolarEclipse_L7390(double x, double y, double igday, int gmonth, int gyear, double tm, double glong, double glat, double hp) + { + var paa = EcRA(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var qaa = EcDec(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var xaa = RightAscensionToHourAngle(DecimalDegreesToDegreeHours(paa), 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var pbb = ParallaxHA(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var qbb = ParallaxDec(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var xbb = HourAngleToRightAscension(pbb, 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var p = (EQELong(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + var q = (EQELat(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + + return (paa, qaa, xaa, pbb, qbb, xbb, p, q); + } + + /// + /// Calculate time of first contact for solar eclipse (UT) + /// + /// + /// Original macro name: UTFirstContactSolarEclipse + /// + public static double UTFirstContactSolarEclipse(double dy, int mn, int yr, int ds, int zc, double glong, double glat) + { + var tp = 2.0 * Math.PI; + + if (SolarEclipseOccurrence(ds, zc, dy, mn, yr).Equals("No solar eclipse")) + return -99.0; + + var dj = NewMoon(ds, zc, dy, mn, yr); + var gday = JulianDateDay(dj); + var gmonth = JulianDateMonth(dj); + var gyear = JulianDateYear(dj); + var igday = gday.Floor(); + var xi = gday - igday; + var utnm = xi * 24.0; + var ut = utnm - 1.0; + var ly = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var my = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var by = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hy = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + ut = utnm + 1.0; + var sb = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians() - ly; + var mz = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var bz = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hz = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + + if (sb < 0.0) + sb = sb + tp; + + var xh = utnm; + var x = my; + var y = by; + var tm = xh - 1.0; + var hp = hy; + var l7390result1 = UTFirstContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + my = l7390result1.p; + by = l7390result1.q; + x = mz; + y = bz; + tm = xh + 1.0; + hp = hz; + var l7390result2 = UTFirstContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + mz = l7390result2.p; + bz = l7390result2.q; + + var x0 = xh + 1.0 - (2.0 * bz / (bz - by)); + var dm = mz - my; + + if (dm < 0.0) + dm = dm + tp; + + var lj = (dm - sb) / 2.0; + var mr = my + (dm * (x0 - xh + 1.0) / 2.0); + ut = x0 - 0.13851852; + var rr = SunDist(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear); + var sr = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + sr = sr + (NutatLong(igday, gmonth, gyear) - 0.00569).ToRadians(); + x = sr; + y = 0.0; + tm = ut; + hp = 0.00004263452 / rr; + var l7390result3 = UTFirstContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + sr = l7390result3.p; + by = by - l7390result3.q; + bz = bz - l7390result3.q; + var p3 = 0.00004263; + var zh = (sr - mr) / lj; + var tc = x0 + zh; + var sh = (((bz - by) * (tc - xh - 1.0) / 2.0) + bz) / lj; + var s2 = sh * sh; + var z2 = zh * zh; + var ps = p3 / (rr * lj); + var z1 = (zh * z2 / (z2 + s2)) + x0; + var h0 = (hy + hz) / (2.0 * lj); + var rm = 0.272446 * h0; + var rn = 0.00465242 / (lj * rr); + var hd = h0 * 0.99834; + var _ru = (hd - rn + ps) * 1.02; + var _rp = (hd + rn + ps) * 1.02; + var pj = Math.Abs(sh * zh / (s2 + z2).SquareRoot()); + var r = rm + rn; + var dd = z1 - x0; + dd = dd * dd - ((z2 - (r * r)) * dd / zh); + + if (dd < 0.0) + return -99.0; + + var zd = dd.SquareRoot(); + var z6 = z1 - zd; + + if (z6 < 0.0) + z6 = z6 + 24.0; + + return z6; + } + + /// + /// Helper function for UTFirstContactSolarEclipse + /// + public static (double paa, double qaa, double xaa, double pbb, double qbb, double xbb, double p, double q) UTFirstContactSolarEclipse_L7390(double x, double y, double igday, int gmonth, int gyear, double tm, double glong, double glat, double hp) + { + var paa = EcRA(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var qaa = EcDec(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var xaa = RightAscensionToHourAngle(DecimalDegreesToDegreeHours(paa), 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var pbb = ParallaxHA(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var qbb = ParallaxDec(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var xbb = HourAngleToRightAscension(pbb, 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var p = (EQELong(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + var q = (EQELat(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + + return (paa, qaa, xaa, pbb, qbb, xbb, p, q); + } + + /// + /// Calculate time of last contact for solar eclipse (UT) + /// + /// + /// Original macro name: UTLastContactSolarEclipse + /// + public static double UTLastContactSolarEclipse(double dy, int mn, int yr, int ds, int zc, double glong, double glat) + { + var tp = 2.0 * Math.PI; + + if (SolarEclipseOccurrence(ds, zc, dy, mn, yr).Equals("No solar eclipse")) + return -99.0; + + var dj = NewMoon(ds, zc, dy, mn, yr); + var gday = JulianDateDay(dj); + var gmonth = JulianDateMonth(dj); + var gyear = JulianDateYear(dj); + var igday = gday.Floor(); + var xi = gday - igday; + var utnm = xi * 24.0; + var ut = utnm - 1.0; + var ly = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var my = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var by = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hy = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + ut = utnm + 1.0; + var sb = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians() - ly; + var mz = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var bz = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hz = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + + if (sb < 0.0) + sb = sb + tp; + + var xh = utnm; + var x = my; + var y = by; + var tm = xh - 1.0; + var hp = hy; + var l7390result1 = UTLastContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + my = l7390result1.p; + by = l7390result1.q; + x = mz; + y = bz; + tm = xh + 1.0; + hp = hz; + var l7390result2 = UTLastContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + mz = l7390result2.p; + bz = l7390result2.q; + + var x0 = xh + 1.0 - (2.0 * bz / (bz - by)); + var dm = mz - my; + + if (dm < 0.0) + dm = dm + tp; + + var lj = (dm - sb) / 2.0; + var mr = my + (dm * (x0 - xh + 1.0) / 2.0); + ut = x0 - 0.13851852; + var rr = SunDist(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear); + var sr = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + sr = sr + (NutatLong(igday, gmonth, gyear) - 0.00569).ToRadians(); + x = sr; + y = 0.0; + tm = ut; + hp = 0.00004263452 / rr; + var l7390result3 = UTLastContactSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + sr = l7390result3.p; + by = by - l7390result3.q; + bz = bz - l7390result3.q; + var p3 = 0.00004263; + var zh = (sr - mr) / lj; + var tc = x0 + zh; + var sh = (((bz - by) * (tc - xh - 1.0) / 2.0) + bz) / lj; + var s2 = sh * sh; + var z2 = zh * zh; + var ps = p3 / (rr * lj); + var z1 = (zh * z2 / (z2 + s2)) + x0; + var h0 = (hy + hz) / (2.0 * lj); + var rm = 0.272446 * h0; + var rn = 0.00465242 / (lj * rr); + var hd = h0 * 0.99834; + var _ru = (hd - rn + ps) * 1.02; + var _rp = (hd + rn + ps) * 1.02; + var pj = Math.Abs(sh * zh / (s2 + z2).SquareRoot()); + var r = rm + rn; + var dd = z1 - x0; + dd = dd * dd - ((z2 - (r * r)) * dd / zh); + + if (dd < 0.0) + return -99.0; + + var zd = dd.SquareRoot(); + var z7 = z1 + zd - Lint((z1 + zd) / 24.0) * 24.0; + + return z7; + } + + /// + /// Helper function for ut_last_contact_solar_eclipse + /// + public static (double paa, double qaa, double xaa, double pbb, double qbb, double xbb, double p, double q) UTLastContactSolarEclipse_L7390(double x, double y, double igday, int gmonth, int gyear, double tm, double glong, double glat, double hp) + { + var paa = EcRA(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var qaa = EcDec(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var xaa = RightAscensionToHourAngle(DecimalDegreesToDegreeHours(paa), 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var pbb = ParallaxHA(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var qbb = ParallaxDec(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var xbb = HourAngleToRightAscension(pbb, 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var p = (EQELong(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + var q = (EQELat(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + + return (paa, qaa, xaa, pbb, qbb, xbb, p, q); + } + + /// + /// Calculate magnitude of solar eclipse. + /// + /// + /// Original macro name: MagSolarEclipse + /// + public static double MagSolarEclipse(double dy, int mn, int yr, int ds, int zc, double glong, double glat) + { + var tp = 2.0 * Math.PI; + + if (SolarEclipseOccurrence(ds, zc, dy, mn, yr).Equals("No solar eclipse")) + return -99.0; + + var dj = NewMoon(ds, zc, dy, mn, yr); + var gday = JulianDateDay(dj); + var gmonth = JulianDateMonth(dj); + var gyear = JulianDateYear(dj); + var igday = gday.Floor(); + var xi = gday - igday; + var utnm = xi * 24.0; + var ut = utnm - 1.0; + var ly = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var my = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var by = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hy = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + ut = utnm + 1.0; + var sb = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians() - ly; + var mz = (MoonLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var bz = (MoonLat(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + var hz = (MoonHP(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + + if (sb < 0.0) + sb = sb + tp; + + var xh = utnm; + var x = my; + var y = by; + var tm = xh - 1.0; + var hp = hy; + var l7390result1 = MagSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + my = l7390result1.p; + by = l7390result1.q; + x = mz; + y = bz; + tm = xh + 1.0; + hp = hz; + var l7390result2 = MagSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + mz = l7390result2.p; + bz = l7390result2.q; + + var x0 = xh + 1.0 - (2.0 * bz / (bz - by)); + var dm = mz - my; + + if (dm < 0.0) + dm = dm + tp; + + var lj = (dm - sb) / 2.0; + var mr = my + (dm * (x0 - xh + 1.0) / 2.0); + ut = x0 - 0.13851852; + var rr = SunDist(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear); + var sr = (SunLong(ut, 0.0, 0.0, 0, 0, igday, gmonth, gyear)).ToRadians(); + sr = sr + (NutatLong(igday, gmonth, gyear) - 0.00569).ToRadians(); + x = sr; + y = 0.0; + tm = ut; + hp = 0.00004263452 / rr; + var l7390result3 = MagSolarEclipse_L7390(x, y, igday, gmonth, gyear, tm, glong, glat, hp); + sr = l7390result3.p; + by = by - l7390result3.q; + bz = bz - l7390result3.q; + var p3 = 0.00004263; + var zh = (sr - mr) / lj; + var tc = x0 + zh; + var sh = (((bz - by) * (tc - xh - 1.0) / 2.0) + bz) / lj; + var s2 = sh * sh; + var z2 = zh * zh; + var ps = p3 / (rr * lj); + var z1 = (zh * z2 / (z2 + s2)) + x0; + var h0 = (hy + hz) / (2.0 * lj); + var rm = 0.272446 * h0; + var rn = 0.00465242 / (lj * rr); + var hd = h0 * 0.99834; + var _ru = (hd - rn + ps) * 1.02; + var _rp = (hd + rn + ps) * 1.02; + var pj = Math.Abs(sh * zh / (s2 + z2).SquareRoot()); + var r = rm + rn; + var dd = z1 - x0; + dd = dd * dd - ((z2 - (r * r)) * dd / zh); + + if (dd < 0.0) + return -99.0; + + var zd = dd.SquareRoot(); + + var mg = (rm + rn - pj) / (2.0 * rn); + + return mg; + } + + /// + /// Helper function for mag_solar_eclipse + /// + public static (double paa, double qaa, double xaa, double pbb, double qbb, double xbb, double p, double q) MagSolarEclipse_L7390(double x, double y, double igday, int gmonth, int gyear, double tm, double glong, double glat, double hp) + { + var paa = EcRA(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var qaa = EcDec(Degrees(x), 0.0, 0.0, Degrees(y), 0.0, 0.0, igday, gmonth, gyear); + var xaa = RightAscensionToHourAngle(DecimalDegreesToDegreeHours(paa), 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var pbb = ParallaxHA(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var qbb = ParallaxDec(xaa, 0.0, 0.0, qaa, 0.0, 0.0, "true", glat, 0.0, Degrees(hp)); + var xbb = HourAngleToRightAscension(pbb, 0.0, 0.0, tm, 0.0, 0.0, 0, 0, igday, gmonth, gyear, glong); + var p = (EQELong(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + var q = (EQELat(xbb, 0.0, 0.0, qbb, 0.0, 0.0, igday, gmonth, gyear)).ToRadians(); + + return (paa, qaa, xaa, pbb, qbb, xbb, p, q); + } } } \ No newline at end of file diff --git a/README.md b/README.md index 131909f..aed7147 100644 --- a/README.md +++ b/README.md @@ -65,4 +65,4 @@ If you're interested in this topic, please buy the book! It provides far more de ### Eclipses - [x] Calculate -> Lunar eclipse occurrence and circumstances -- [ ] Calculate -> Solar eclipse occurrence and circumstances +- [x] Calculate -> Solar eclipse occurrence and circumstances