Merge branch 'main' into feature/dusk/ConvertEclipticCoordinateToEqua…

…torial
observerly · Feb 21, 2022 · 6e3eee2 · 6e3eee2
2 parents f048a66 + 4af4981
commit 6e3eee2
Show file tree

Hide file tree

Showing 2 changed files with 259 additions and 0 deletions.
diff --git a/pkg/dusk/lunar.go b/pkg/dusk/lunar.go
@@ -196,6 +196,242 @@ func GetLunarEquatorialPosition(datetime time.Time) EquatorialCoordinate {
 	}
 }
 
+/*
+	GetMoonEclipticPosition()
+
+	N.B. the ecliptic position is referenced to mean equinox of date and do not include the effect of nutation.
+
+	@param datetime - the datetime in UTC of the observer
+	@returns the geocentric ecliptic coodinate (λ - geocentric longitude, β - geocentric latidude and Δ distance between centers of the Earth and Moon, in km) of the Moon.
+*/
+func GetLunarEclipticPosition(datetime time.Time) EclipticCoordinate {
+	var T float64 = GetCurrentJulianCenturyRelativeToJ2000(datetime)
+
+	var D float64 = GetLunarMeanElongation(T)
+
+	var Lʹ float64 = GetLunarMeanLongitude(T)
+
+	var M float64 = GetSolarMeanAnomaly(T)
+
+	var Mʹ float64 = GetLunarMeanAnomaly(T)
+
+	var F float64 = GetLunarArgumentOfLatitude(T)
+
+	var A1 float64 = 119.75 + 131.849*T
+
+	var A2 float64 = 53.09 + 479264.29*T
+
+	var A3 float64 = 313.45 + 481266.484*T
+
+	var E float64 = 1 - 0.002516*T - 0.0000074*math.Pow(T, 2)
+
+	var E2 float64 = math.Pow(E, 2)
+
+	var Σl float64 = 3958*sinx(A1) + 1962*sinx(Lʹ-F) + 318*sinx(A2)
+
+	var Σr float64 = 0
+
+	var Σb float64 = -2235*sinx(Lʹ) + 382*sinx(A3) + 175*sinx(A1-F) + 175*sinx(A1+F) + 127*sinx(Lʹ-Mʹ) - 115*sinx(Lʹ+Mʹ)
+
+	for i := range ta {
+		r := &ta[i]
+
+		sa, ca := sincosx(D*r.D + M*r.M + Mʹ*r.Mʹ + F*r.F)
+
+		switch r.M {
+		case 0:
+			Σl += r.Σl * sa
+			Σr += r.Σr * ca
+		case 1, -1:
+			Σl += r.Σl * sa * E
+			Σr += r.Σr * ca * E
+		case 2, -2:
+			Σl += r.Σl * sa * E2
+			Σr += r.Σr * ca * E2
+		}
+	}
+
+	for i := range tb {
+		r := &tb[i]
+
+		sb := sinx(D*r.D + M*r.M + Mʹ*r.Mʹ + F*r.F)
+
+		switch r.M {
+		case 0:
+			Σb += r.Σb * sb
+		case 1, -1:
+			Σb += r.Σb * sb * E
+		case 2, -2:
+			Σb += r.Σb * sb * E2
+		}
+	}
+
+	return EclipticCoordinate{
+		λ: Lʹ + Σl/1000000,
+		β: Σb / 1000000,
+		Δ: 385000.56 + Σr/1000,
+	}
+}
+
+type tas struct{ D, M, Mʹ, F, Σl, Σr float64 }
+
+var ta = [...]tas{
+	{0, 0, 1, 0, 6288774, -20905355},
+	{2, 0, -1, 0, 1274027, -3699111},
+	{2, 0, 0, 0, 658314, -2955968},
+	{0, 0, 2, 0, 213618, -569925},
+
+	{0, 1, 0, 0, -185116, 48888},
+	{0, 0, 0, 2, -114332, -3149},
+	{2, 0, -2, 0, 58793, 246158},
+	{2, -1, -1, 0, 57066, -152138},
+
+	{2, 0, 1, 0, 53322, -170733},
+	{2, -1, 0, 0, 45758, -204586},
+	{0, 1, -1, 0, -40923, -129620},
+	{1, 0, 0, 0, -34720, 108743},
+
+	{0, 1, 1, 0, -30383, 104755},
+	{2, 0, 0, -2, 15327, 10321},
+	{0, 0, 1, 2, -12528, 0},
+	{0, 0, 1, -2, 10980, 79661},
+
+	{4, 0, -1, 0, 10675, -34782},
+	{0, 0, 3, 0, 10034, -23210},
+	{4, 0, -2, 0, 8548, -21636},
+	{2, 1, -1, 0, -7888, 24208},
+
+	{2, 1, 0, 0, -6766, 30824},
+	{1, 0, -1, 0, -5163, -8379},
+	{1, 1, 0, 0, 4987, -16675},
+	{2, -1, 1, 0, 4036, -12831},
+
+	{2, 0, 2, 0, 3994, -10445},
+	{4, 0, 0, 0, 3861, -11650},
+	{2, 0, -3, 0, 3665, 14403},
+	{0, 1, -2, 0, -2689, -7003},
+
+	{2, 0, -1, 2, -2602, 0},
+	{2, -1, -2, 0, 2390, 10056},
+	{1, 0, 1, 0, -2348, 6322},
+	{2, -2, 0, 0, 2236, -9884},
+
+	{0, 1, 2, 0, -2120, 5751},
+	{0, 2, 0, 0, -2069, 0},
+	{2, -2, -1, 0, 2048, -4950},
+	{2, 0, 1, -2, -1773, 4130},
+
+	{2, 0, 0, 2, -1595, 0},
+	{4, -1, -1, 0, 1215, -3958},
+	{0, 0, 2, 2, -1110, 0},
+	{3, 0, -1, 0, -892, 3258},
+
+	{2, 1, 1, 0, -810, 2616},
+	{4, -1, -2, 0, 759, -1897},
+	{0, 2, -1, 0, -713, -2117},
+	{2, 2, -1, 0, -700, 2354},
+
+	{2, 1, -2, 0, 691, 0},
+	{2, -1, 0, -2, 596, 0},
+	{4, 0, 1, 0, 549, -1423},
+	{0, 0, 4, 0, 537, -1117},
+
+	{4, -1, 0, 0, 520, -1571},
+	{1, 0, -2, 0, -487, -1739},
+	{2, 1, 0, -2, -399, 0},
+	{0, 0, 2, -2, -381, -4421},
+
+	{1, 1, 1, 0, 351, 0},
+	{3, 0, -2, 0, -340, 0},
+	{4, 0, -3, 0, 330, 0},
+	{2, -1, 2, 0, 327, 0},
+
+	{0, 2, 1, 0, -323, 1165},
+	{1, 1, -1, 0, 299, 0},
+	{2, 0, 3, 0, 294, 0},
+	{2, 0, -1, -2, 0, 8752},
+}
+
+type tbs struct{ D, M, Mʹ, F, Σb float64 }
+
+var tb = [...]tbs{
+	{0, 0, 0, 1, 5128122},
+	{0, 0, 1, 1, 280602},
+	{0, 0, 1, -1, 277693},
+	{2, 0, 0, -1, 173237},
+
+	{2, 0, -1, 1, 55413},
+	{2, 0, -1, -1, 46271},
+	{2, 0, 0, 1, 32573},
+	{0, 0, 2, 1, 17198},
+
+	{2, 0, 1, -1, 9266},
+	{0, 0, 2, -1, 8822},
+	{2, -1, 0, -1, 8216},
+	{2, 0, -2, -1, 4324},
+
+	{2, 0, 1, 1, 4200},
+	{2, 1, 0, -1, -3359},
+	{2, -1, -1, 1, 2463},
+	{2, -1, 0, 1, 2211},
+
+	{2, -1, -1, -1, 2065},
+	{0, 1, -1, -1, -1870},
+	{4, 0, -1, -1, 1828},
+	{0, 1, 0, 1, -1794},
+
+	{0, 0, 0, 3, -1749},
+	{0, 1, -1, 1, -1565},
+	{1, 0, 0, 1, -1491},
+	{0, 1, 1, 1, -1475},
+
+	{0, 1, 1, -1, -1410},
+	{0, 1, 0, -1, -1344},
+	{1, 0, 0, -1, -1335},
+	{0, 0, 3, 1, 1107},
+
+	{4, 0, 0, -1, 1021},
+	{4, 0, -1, 1, 833},
+
+	{0, 0, 1, -3, 777},
+	{4, 0, -2, 1, 671},
+	{2, 0, 0, -3, 607},
+	{2, 0, 2, -1, 596},
+
+	{2, -1, 1, -1, 491},
+	{2, 0, -2, 1, -451},
+	{0, 0, 3, -1, 439},
+	{2, 0, 2, 1, 422},
+
+	{2, 0, -3, -1, 421},
+	{2, 1, -1, 1, -366},
+	{2, 1, 0, 1, -351},
+	{4, 0, 0, 1, 331},
+
+	{2, -1, 1, 1, 315},
+	{2, -2, 0, -1, 302},
+	{0, 0, 1, 3, -283},
+	{2, 1, 1, -1, -229},
+
+	{1, 1, 0, -1, 223},
+	{1, 1, 0, 1, 223},
+	{0, 1, -2, -1, -220},
+	{2, 1, -1, -1, -220},
+
+	{1, 0, 1, 1, -185},
+	{2, -1, -2, -1, 181},
+	{0, 1, 2, 1, -177},
+	{4, 0, -2, -1, 176},
+
+	{4, -1, -1, -1, 166},
+	{1, 0, 1, -1, -164},
+	{4, 0, 1, -1, 132},
+	{1, 0, -1, -1, -119},
+
+	{4, -1, 0, -1, 115},
+	{2, -2, 0, 1, 107},
+}
+
 /*
 	GetLunarHourAngle()
 

diff --git a/pkg/dusk/lunar_test.go b/pkg/dusk/lunar_test.go
@@ -172,6 +172,29 @@ func TestGetLunarEquatorialPositionDeclination(t *testing.T) {
 		t.Errorf("got %f, wanted %f", got, want)
 	}
 }
+func TestGetLunarEclipticPositionLongitude(t *testing.T) {
+	var ec EclipticCoordinate = GetLunarEclipticPosition(d)
+
+	var got float64 = ec.λ
+
+	var want float64 = 133.162655
+
+	if math.Abs(got-want) > 0.15 {
+		t.Errorf("got %f, wanted %f", got, want)
+	}
+}
+
+func TestGetLunarEclipticPositionLatitude(t *testing.T) {
+	var ec EclipticCoordinate = GetLunarEclipticPosition(d)
+
+	var got float64 = ec.β
+
+	var want float64 = -3.229126
+
+	if math.Abs(got-want) > 0.1 {
+		t.Errorf("got %f, wanted %f", got, want)
+	}
+}
 
 func TestGetLunarHourAngle(t *testing.T) {
 	var eq EquatorialCoordinate = GetLunarEquatorialPosition(datetime)