// Nutation returns nutation in longitude (Δψ) and nutation in obliquity (Δε)
// for a given JDE.
//
// JDE = UT + ΔT, see package deltat.
//
// Computation is by 1980 IAU theory, with terms < .0003″ neglected.
func Nutation(jde float64) (Δψ, Δε unit.Angle) {
	T := base.J2000Century(jde)
	D := base.Horner(T,
		297.85036, 445267.11148, -0.0019142, 1./189474) * math.Pi / 180
	M := base.Horner(T,
		357.52772, 35999.050340, -0.0001603, -1./300000) * math.Pi / 180
	N := base.Horner(T,
		134.96298, 477198.867398, 0.0086972, 1./5620) * math.Pi / 180
	F := base.Horner(T,
		93.27191, 483202.017538, -0.0036825, 1./327270) * math.Pi / 180
	Ω := base.Horner(T,
		125.04452, -1934.136261, 0.0020708, 1./450000) * math.Pi / 180
	// sum in reverse order to accumulate smaller terms first
	var Δψs, Δεs float64
	for i := len(table22A) - 1; i >= 0; i-- {
		row := table22A[i]
		arg := row.d*D + row.m*M + row.n*N + row.f*F + row.ω*Ω
		s, c := math.Sincos(arg)
		Δψs += s * (row.s0 + row.s1*T)
		Δεs += c * (row.c0 + row.c1*T)
	}
	Δψ = unit.AngleFromSec(Δψs * .0001)
	Δε = unit.AngleFromSec(Δεs * .0001)
	return
}

// TrueVSOP87 returns the true geometric position of the sun as ecliptic coordinates.
//
// Result computed by full VSOP87 theory.  Result is at equator and equinox
// of date in the FK5 frame.  It does not include nutation or aberration.
//
//	s: ecliptic longitude
//	β: ecliptic latitude
//	R: range in AU
func TrueVSOP87(e *pp.V87Planet, jde float64) (s, β unit.Angle, R float64) {
	l, b, r := e.Position(jde)
	s = l + math.Pi
	// FK5 correction.
	λp := base.Horner(base.J2000Century(jde),
		s.Rad(), -1.397*math.Pi/180, -.00031*math.Pi/180)
	sλp, cλp := math.Sincos(λp)
	Δβ := unit.AngleFromSec(.03916).Mul(cλp - sλp)
	// (25.9) p. 166
	s -= unit.AngleFromSec(.09033)
	return s.Mod1(), Δβ - b, r
}

// ApproxNutation returns a fast approximation of nutation in longitude (Δψ)
// and nutation in obliquity (Δε) for a given JDE.
//
// Accuracy is 0.5″ in Δψ, 0.1″ in Δε.
func ApproxNutation(jde float64) (Δψ, Δε unit.Angle) {
	T := (jde - base.J2000) / 36525
	Ω := (125.04452 - 1934.136261*T) * math.Pi / 180
	L := (280.4665 + 36000.7698*T) * math.Pi / 180
	N := (218.3165 + 481267.8813*T) * math.Pi / 180
	sΩ, cΩ := math.Sincos(Ω)
	s2L, c2L := math.Sincos(2 * L)
	s2N, c2N := math.Sincos(2 * N)
	s2Ω, c2Ω := math.Sincos(2 * Ω)
	Δψ = unit.AngleFromSec(-17.2*sΩ - 1.32*s2L - 0.23*s2N + 0.21*s2Ω)
	Δε = unit.AngleFromSec(9.2*cΩ + 0.57*c2L + 0.1*c2N - 0.09*c2Ω)
	return
}

// mn as separate function for testing purposes
func mn(epochFrom, epochTo float64) (m, nα unit.HourAngle, nδ unit.Angle) {
	T := (epochTo - epochFrom) * .01
	m = unit.HourAngleFromSec(3.07496 + 0.00186*T)
	nα = unit.HourAngleFromSec(1.33621 - 0.00057*T)
	nδ = unit.AngleFromSec(20.0431 - 0.0085*T)
	return
}

func ExampleProperMotion3D() {
	// Example 21.d, p. 141.
	eqFrom := &coord.Equatorial{
		RA:  unit.NewRA(6, 45, 8.871),
		Dec: unit.NewAngle('-', 16, 42, 57.99),
	}
	mra := unit.HourAngleFromSec(-0.03847)
	mdec := unit.AngleFromSec(-1.2053)
	r := 2.64           // given in correct unit
	mr := -7.6 / 977792 // magic conversion factor
	eqTo := &coord.Equatorial{}
	fmt.Printf("Δr = %.9f, Δα = %.10f, Δδ = %.10f\n", mr, mra, mdec)
	for _, epoch := range []float64{1000, 0, -1000, -2000, -10000} {
		precess.ProperMotion3D(eqFrom, eqTo, 2000, epoch, r, mr, mra, mdec)
		fmt.Printf("%8.1f  %0.2d  %0.1d\n", epoch,
			sexa.FmtRA(eqTo.RA), sexa.FmtAngle(eqTo.Dec))
	}
	// Output:
	// Δr = -0.000007773, Δα = -0.0000027976, Δδ = -0.0000058435
	//   1000.0  6ʰ45ᵐ47ˢ.16  -16°22′56″.0
	//      0.0  6ʰ46ᵐ25ˢ.09  -16°03′00″.8
	//  -1000.0  6ʰ47ᵐ02ˢ.67  -15°43′12″.3
	//  -2000.0  6ʰ47ᵐ39ˢ.91  -15°23′30″.6
	// -10000.0  6ʰ52ᵐ25ˢ.72  -12°50′06″.7
}

// Exercise, p. 136.
func TestPosition(t *testing.T) {
	eqFrom := &coord.Equatorial{
		unit.NewRA(2, 31, 48.704),
		unit.NewAngle(' ', 89, 15, 50.72),
	}
	eqTo := &coord.Equatorial{}
	mα := unit.HourAngleFromSec(0.19877)
	mδ := unit.AngleFromSec(-0.0152)
	for _, tc := range []struct {
		α, δ string
		jde  float64
	}{
		{"1ʰ22ᵐ33.90ˢ", "88°46′26.18″", base.BesselianYearToJDE(1900)},
		{"3ʰ48ᵐ16.43ˢ", "89°27′15.38″", base.JulianYearToJDE(2050)},
		{"5ʰ53ᵐ29.17ˢ", "89°32′22.18″", base.JulianYearToJDE(2100)},
	} {
		epochTo := base.JDEToJulianYear(tc.jde)
		precess.Position(eqFrom, eqTo, 2000.0, epochTo, mα, mδ)
		αStr := fmt.Sprintf("%.2s", sexa.FmtRA(eqTo.RA))
		δStr := fmt.Sprintf("%.2s", sexa.FmtAngle(eqTo.Dec))
		if αStr != tc.α {
			t.Fatal("got:", αStr, "want:", tc.α)
		}
		if δStr != tc.δ {
			t.Fatal(δStr)
		}
	}
}

func TestPrecessor_Precess(t *testing.T) {
	// Exercise, p. 136.
	eqFrom := &coord.Equatorial{
		RA:  unit.NewRA(2, 31, 48.704),
		Dec: unit.NewAngle(' ', 89, 15, 50.72),
	}
	mα := unit.HourAngleFromSec(.19877)
	mδ := unit.AngleFromSec(-.0152)
	epochs := []float64{
		base.JDEToJulianYear(base.B1900),
		2050,
		2100,
	}
	answer := []string{
		"α = 1ʰ22ᵐ33ˢ.90   δ = +88°46′26″.18",
		"α = 3ʰ48ᵐ16ˢ.43   δ = +89°27′15″.38",
		"α = 5ʰ53ᵐ29ˢ.17   δ = +89°32′22″.18",
	}
	eqTo := &coord.Equatorial{}
	for i, epochTo := range epochs {
		precess.Position(eqFrom, eqTo, 2000, epochTo, mα, mδ)
		if answer[i] != fmt.Sprintf("α = %0.2d   δ = %+0.2d",
			sexa.FmtRA(eqTo.RA), sexa.FmtAngle(eqTo.Dec)) {
			t.Fatal(i)
		}
	}
}

// MeanObliquity returns mean obliquity (ε₀) following the IAU 1980
// polynomial.
//
// Accuracy is 1″ over the range 1000 to 3000 years and 10″ over the range
// 0 to 4000 years.
func MeanObliquity(jde float64) unit.Angle {
	// (22.2) p. 147
	return unit.AngleFromSec(base.Horner(base.J2000Century(jde),
		unit.FromSexaSec(' ', 23, 26, 21.448),
		-46.815,
		-0.00059,
		0.001813))
}

// "rectangular coordinate" solution, p. 113.
func TestMinSepRect(t *testing.T) {
	sep, err := angle.MinSepRect(jd1, jd3, r1, d1, r2, d2)
	if err != nil {
		t.Fatal(err)
	}
	answer := unit.AngleFromSec(224) // on p. 111
	if math.Abs((sep-answer).Rad()/sep.Rad()) > 1e-2 {
		t.Fatal(sep, answer)
	}

}

// Test with proper motion of Regulus, with equatorial motions given
// in Example 21.a, p. 132, and ecliptic motions given in table 21.A,
// p. 138.
func TestEqProperMotionToEcl(t *testing.T) {
	ε := coord.NewObliquity(nutation.MeanObliquity(base.J2000))
	mλ, mβ := eqProperMotionToEcl(
		// eq motions from p. 132.
		unit.NewHourAngle('-', 0, 0, 0.0169),
		unit.NewAngle(' ', 0, 0, 0.006),
		2000.0,
		// eq coordinates from p. 132.
		new(coord.Ecliptic).EqToEcl(&coord.Equatorial{
			RA:  unit.NewRA(10, 8, 22.3),
			Dec: unit.NewAngle(' ', 11, 58, 2),
		}, ε))
	d := math.Abs((mλ - unit.AngleFromSec(-.2348)).Rad() / mλ.Rad())
	if d*169 > 1 { // 169 = significant digits of given lon
		t.Fatal("mλ")
	}
	d = math.Abs((mβ - unit.AngleFromSec(-.0813)).Rad() / mβ.Rad())
	if d*6 > 1 { // 6 = significant digit of given lat
		t.Fatal("mβ")
	}
}

// perigee parallax
func (l *la) pp() unit.Angle {
	return unit.AngleFromSec(
		3629.215 +
			63.224*math.Cos(2*l.D) +
			-6.990*math.Cos(4*l.D) +
			(2.834-0.0071*l.T)*math.Cos(2*l.D-l.M) +
			1.927*math.Cos(6*l.D) +
			-1.263*math.Cos(l.D) +
			-0.702*math.Cos(8*l.D) +
			(0.696-0.0017*l.T)*math.Cos(l.M) +
			-0.690*math.Cos(2*l.F) +
			(-0.629+0.0016*l.T)*math.Cos(4*l.D-l.M) +
			-0.392*math.Cos(2*(l.D-l.F)) +
			0.297*math.Cos(10*l.D) +
			0.260*math.Cos(6*l.D-l.M) +
			0.201*math.Cos(3*l.D) +
			-0.161*math.Cos(2*l.D+l.M) +
			0.157*math.Cos(l.D+l.M) +
			-0.138*math.Cos(12*l.D) +
			-0.127*math.Cos(8*l.D-l.M) +
			0.104*math.Cos(2*(l.D+l.F)) +
			0.104*math.Cos(2*(l.D-l.M)) +
			-0.079*math.Cos(5*l.D) +
			0.068*math.Cos(14*l.D) +
			0.067*math.Cos(10*l.D-l.M) +
			0.054*math.Cos(4*l.D+l.M) +
			-0.038*math.Cos(12*l.D-l.M) +
			-0.038*math.Cos(4*l.D-2*l.M) +
			0.037*math.Cos(7*l.D) +
			-0.037*math.Cos(4*l.D+2*l.F) +
			-0.035*math.Cos(16*l.D) +
			-0.030*math.Cos(3*l.D+l.M) +
			0.029*math.Cos(l.D-l.M) +
			-0.025*math.Cos(6*l.D+l.M) +
			0.023*math.Cos(2*l.M) +
			0.023*math.Cos(14*l.D-l.M) +
			-0.023*math.Cos(2*(l.D+l.M)) +
			0.022*math.Cos(6*l.D-2*l.M) +
			-0.021*math.Cos(2*l.D-2*l.F-l.M) +
			-0.020*math.Cos(9*l.D) +
			0.019*math.Cos(18*l.D) +
			0.017*math.Cos(6*l.D+2*l.F) +
			0.014*math.Cos(2*l.F-l.M) +
			-0.014*math.Cos(16*l.D-l.M) +
			0.013*math.Cos(4*l.D-2*l.F) +
			0.012*math.Cos(8*l.D+l.M) +
			0.011*math.Cos(11*l.D) +
			0.010*math.Cos(5*l.D+l.M) +
			-0.010*math.Cos(20*l.D))
}

// MeanObliquityLaskar returns mean obliquity (ε₀) following the Laskar
// 1986 polynomial.
//
// Accuracy over the range 1000 to 3000 years is .01″.
//
// Accuracy over the valid date range of -8000 to +12000 years is
// "a few seconds."
func MeanObliquityLaskar(jde float64) unit.Angle {
	// (22.3) p. 147
	return unit.AngleFromSec(base.Horner(base.J2000Century(jde)*.01,
		unit.FromSexaSec(' ', 23, 26, 21.448),
		-4680.93,
		-1.55,
		1999.25,
		-51.38,
		-249.67,
		-39.05,
		7.12,
		27.87,
		5.79,
		2.45))
}

func ExampleStellar() {
	// Exercise, p. 119.
	day1 := 7.
	day5 := 27.
	r2 := []unit.Angle{
		unit.NewRA(15, 3, 51.937).Angle(),
		unit.NewRA(15, 9, 57.327).Angle(),
		unit.NewRA(15, 15, 37.898).Angle(),
		unit.NewRA(15, 20, 50.632).Angle(),
		unit.NewRA(15, 25, 32.695).Angle(),
	}
	d2 := []unit.Angle{
		unit.NewAngle('-', 8, 57, 34.51),
		unit.NewAngle('-', 9, 9, 03.88),
		unit.NewAngle('-', 9, 17, 37.94),
		unit.NewAngle('-', 9, 23, 16.25),
		unit.NewAngle('-', 9, 26, 01.01),
	}
	jd := julian.CalendarGregorianToJD(1996, 2, 17)
	dt := jd - base.J2000
	dy := dt / base.JulianYear
	dc := dy / 100
	fmt.Printf("%.2f years\n", dy)
	fmt.Printf("%.4f century\n", dc)

	pmr := -.649 // sec/cen
	pmd := -1.91 // sec/cen
	r1 := unit.NewRA(15, 17, 0.421) + unit.RAFromSec(pmr*dc)
	d1 := unit.NewAngle('-', 9, 22, 58.54) + unit.AngleFromSec(pmd*dc)
	fmt.Printf("α′ = %.3d, δ′ = %.2d\n", sexa.FmtRA(r1), sexa.FmtAngle(d1))

	day, dd, err := conjunction.Stellar(day1, day5, r1.Angle(), d1, r2, d2)
	if err != nil {
		fmt.Println(err)
		return
	}
	fmt.Println(sexa.FmtAngle(dd))
	dInt, dFrac := math.Modf(day)
	fmt.Printf("1996 February %d at %s TD\n", int(dInt),
		sexa.FmtTime(unit.TimeFromDay(dFrac)))

	// Output:
	// -3.87 years
	// -0.0387 century
	// α′ = 15ʰ17ᵐ0ˢ.446, δ′ = -9°22′58″.47
	// 3′38″
	// 1996 February 18 at 6ʰ36ᵐ55ˢ TD
}

// E computes the "equation of time" for the given JDE.
//
// Parameter e must be a planetposition.V87Planet object for Earth obtained
// with planetposition.LoadPlanet.
//
// Result is equation of time as an hour angle.
func E(jde float64, e *pp.V87Planet) unit.HourAngle {
	τ := base.J2000Century(jde) * .1
	L0 := l0(τ)
	// code duplicated from solar.ApparentEquatorialVSOP87 so that
	// we can keep Δψ and cε
	s, β, R := solar.TrueVSOP87(e, jde)
	Δψ, Δε := nutation.Nutation(jde)
	a := unit.AngleFromSec(-20.4898).Div(R)
	λ := s + Δψ + a
	ε := nutation.MeanObliquity(jde) + Δε
	sε, cε := ε.Sincos()
	α, _ := coord.EclToEq(λ, β, sε, cε)
	// (28.1) p. 183
	E := L0 - unit.AngleFromDeg(.0057183) - unit.Angle(α) + Δψ.Mul(cε)
	return unit.HourAngle((E + math.Pi).Mod1() - math.Pi)
}

func ExamplePositionRonVondrak() {
	// Example 23.b, p. 156
	jd := julian.CalendarGregorianToJD(2028, 11, 13.19)
	eq := &coord.Equatorial{
		RA:  unit.NewRA(2, 44, 11.986),
		Dec: unit.NewAngle(' ', 49, 13, 42.48),
	}
	apparent.PositionRonVondrak(eq, eq, base.JDEToJulianYear(jd),
		unit.HourAngleFromSec(.03425),
		unit.AngleFromSec(-.0895))
	fmt.Printf("α = %0.3d\n", sexa.FmtRA(eq.RA))
	fmt.Printf("δ = %0.2d\n", sexa.FmtAngle(eq.Dec))
	// Output:
	// α = 2ʰ46ᵐ14ˢ.392
	// δ = 49°21′07″.45
}

func ExampleApproxPosition() {
	// Example 21.a, p. 132.
	eq := &coord.Equatorial{
		unit.NewRA(10, 8, 22.3),
		unit.NewAngle(' ', 11, 58, 2),
	}
	epochFrom := 2000.0
	epochTo := 1978.0
	mα := unit.HourAngleFromSec(-0.0169)
	mδ := unit.AngleFromSec(0.006)
	precess.ApproxPosition(eq, eq, epochFrom, epochTo, mα, mδ)
	fmt.Printf("%0.1d\n", sexa.FmtRA(eq.RA))
	fmt.Printf("%+0d\n", sexa.FmtAngle(eq.Dec))
	// Output:
	// 10ʰ07ᵐ12ˢ.1
	// +12°04′32″
}

func ExamplePosition() {
	// Example 21.b, p. 135.
	eq := &coord.Equatorial{
		unit.NewRA(2, 44, 11.986),
		unit.NewAngle(' ', 49, 13, 42.48),
	}
	epochFrom := 2000.0
	jdTo := julian.CalendarGregorianToJD(2028, 11, 13.19)
	epochTo := base.JDEToJulianYear(jdTo)
	precess.Position(eq, eq, epochFrom, epochTo,
		unit.HourAngleFromSec(0.03425),
		unit.AngleFromSec(-0.0895))
	fmt.Printf("%0.3d\n", sexa.FmtRA(eq.RA))
	fmt.Printf("%+0.2d\n", sexa.FmtAngle(eq.Dec))
	// Output:
	// 2ʰ46ᵐ11ˢ.331
	// +49°20′54″.54
}

func ExamplePosition() {
	// Example 57.a, p. 398
	M := binary.M(1980, 1934.008, 41.623)
	fmt.Printf("M = %.3f\n", M.Deg())
	E, err := kepler.Kepler1(.2763, M, 4)
	if err != nil {
		fmt.Println(err)
		return
	}
	fmt.Printf("E = %.3f\n", E.Deg())
	θ, ρ := binary.Position(.2763, unit.AngleFromSec(.907),
		unit.AngleFromDeg(59.025), unit.AngleFromDeg(23.717),
		unit.AngleFromDeg(219.907), E)
	fmt.Printf("θ = %.1f\n", θ.Deg())
	fmt.Printf("ρ = %.3f\n", ρ.Sec())
	// Output:
	// M = 37.788
	// E = 49.896
	// θ = 318.4
	// ρ = 0.411
}

// apogee parallax
func (l *la) ap() unit.Angle {
	return unit.AngleFromSec(
		3245.251 +
			-9.147*math.Cos(2*l.D) +
			-.841*math.Cos(l.D) +
			.697*math.Cos(2*l.F) +
			(-.656+.0016*l.T)*math.Cos(l.M) +
			.355*math.Cos(4*l.D) +
			.159*math.Cos(2*l.D-l.M) +
			.127*math.Cos(l.D+l.M) +
			.065*math.Cos(4*l.D-l.M) +
			.052*math.Cos(6*l.D) +
			.043*math.Cos(2*l.D+l.M) +
			.031*math.Cos(2*(l.D+l.F)) +
			-.023*math.Cos(2*(l.D-l.F)) +
			.022*math.Cos(2*(l.D-l.M)) +
			.019*math.Cos(2*(l.D+l.M)) +
			-.016*math.Cos(2*l.M) +
			.014*math.Cos(6*l.D-l.M) +
			.01*math.Cos(8*l.D))
}

// Ephemeris returns the apparent orientation of the sun at the given jd.
//
// Results:
//	P:  Position angle of the solar north pole.
//	B0: Heliographic latitude of the center of the solar disk.
//	L0: Heliographic longitude of the center of the solar disk.
func Ephemeris(jd float64, e *pp.V87Planet) (P, B0, L0 unit.Angle) {
	θ := unit.Angle((jd - 2398220) * 2 * math.Pi / 25.38)
	I := unit.AngleFromDeg(7.25)
	K := unit.AngleFromDeg(73.6667) +
		unit.AngleFromDeg(1.3958333).Mul((jd-2396758)/base.JulianCentury)

	L, _, R := solar.TrueVSOP87(e, jd)
	Δψ, Δε := nutation.Nutation(jd)
	ε0 := nutation.MeanObliquity(jd)
	ε := ε0 + Δε
	λ := L - unit.AngleFromSec(20.4898).Div(R)
	λp := λ + Δψ

	sλK, cλK := (λ - K).Sincos()
	sI, cI := I.Sincos()

	tx := -(λp.Cos() * ε.Tan())
	ty := -(cλK * I.Tan())
	P = unit.Angle(math.Atan(tx) + math.Atan(ty))
	B0 = unit.Angle(math.Asin(sλK * sI))
	η := unit.Angle(math.Atan2(-sλK*cI, -cλK))
	L0 = (η - θ).Mod1()
	return
}