// Position returns observed equatorial coordinates of a planet at a given time.
//
// Argument p must be a valid V87Planet object for the observed planet.
// Argument earth must be a valid V87Planet object for Earth.
//
// Results are right ascension and declination, α and δ in radians.
func Position(p, earth *pp.V87Planet, jde float64) (α, δ float64) {
	L0, B0, R0 := earth.Position(jde)
	L, B, R := p.Position(jde)
	sB0, cB0 := math.Sincos(B0)
	sL0, cL0 := math.Sincos(L0)
	sB, cB := math.Sincos(B)
	sL, cL := math.Sincos(L)
	x := R*cB*cL - R0*cB0*cL0
	y := R*cB*sL - R0*cB0*sL0
	z := R*sB - R0*sB0
	{
		Δ := math.Sqrt(x*x + y*y + z*z) // (33.4) p. 224
		τ := base.LightTime(Δ)
		// repeating with jde-τ
		L, B, R = p.Position(jde - τ)
		sB, cB = math.Sincos(B)
		sL, cL = math.Sincos(L)
		x = R*cB*cL - R0*cB0*cL0
		y = R*cB*sL - R0*cB0*sL0
		z = R*sB - R0*sB0
	}
	λ := math.Atan2(y, x)                // (33.1) p. 223
	β := math.Atan2(z, math.Hypot(x, y)) // (33.2) p. 223
	Δλ, Δβ := apparent.EclipticAberration(λ, β, jde)
	λ, β = pp.ToFK5(λ+Δλ, β+Δβ, jde)
	Δψ, Δε := nutation.Nutation(jde)
	λ += Δψ
	sε, cε := math.Sincos(nutation.MeanObliquity(jde) + Δε)
	return coord.EclToEq(λ, β, sε, cε)
	// Meeus gives a formula for elongation but doesn't spell out how to
	// obtaion term λ0 and doesn't give an example solution.
}

// TrueEquatorial returns the true geometric position of the Sun as equatorial coordinates.
func TrueEquatorial(jde float64) (α, δ float64) {
	s, _ := True(base.J2000Century(jde))
	ε := nutation.MeanObliquity(jde)
	ss, cs := math.Sincos(s)
	sε, cε := math.Sincos(ε)
	// (25.6, 25.7) p. 165
	return math.Atan2(cε*ss, cs), sε * ss
}

// ApparentEquatorial returns the apparent position of the Sun as equatorial coordinates.
//
//	α: right ascension in radians
//	δ: declination in radians
func ApparentEquatorial(jde float64) (α, δ float64) {
	T := base.J2000Century(jde)
	λ := ApparentLongitude(T)
	ε := nutation.MeanObliquity(jde)
	sλ, cλ := math.Sincos(λ)
	// (25.8) p. 165
	sε, cε := math.Sincos(ε + .00256*math.Pi/180*math.Cos(node(T)))
	return math.Atan2(cε*sλ, cλ), math.Asin(sε * sλ)
}

func TestIAUvsLaskar(t *testing.T) {
	for _, y := range []int{1000, 2000, 3000} {
		jd := julian.CalendarGregorianToJD(y, 0, 0)
		i := nutation.MeanObliquity(jd)
		l := nutation.MeanObliquityLaskar(jd)
		if math.Abs((i - l).Sec()) > 1 {
			t.Fatal(y)
		}
	}
	for _, y := range []int{0, 4000} {
		jd := julian.CalendarGregorianToJD(y, 0, 0)
		i := nutation.MeanObliquity(jd)
		l := nutation.MeanObliquityLaskar(jd)
		if math.Abs((i - l).Sec()) > 10 {
			t.Fatal(y)
		}
	}
}

// TrueEquatorial returns the true geometric position of the Sun as equatorial coordinates.
func TrueEquatorial(jde float64) (α unit.RA, δ unit.Angle) {
	s, _ := True(base.J2000Century(jde))
	ε := nutation.MeanObliquity(jde)
	ss, cs := s.Sincos()
	sε, cε := ε.Sincos()
	// (25.6, 25.7) p. 165
	α = unit.RAFromRad(math.Atan2(cε*ss, cs))
	δ = unit.Angle(math.Asin(sε * ss))
	return
}

// Position returns rectangular coordinates referenced to the mean equinox of date.
func Position(e *pp.V87Planet, jde float64) (x, y, z float64) {
	// (26.1) p. 171
	s, β, R := solar.TrueVSOP87(e, jde)
	sε, cε := math.Sincos(nutation.MeanObliquity(jde))
	ss, cs := math.Sincos(s)
	sβ := math.Sin(β)
	x = R * cs
	y = R * (ss*cε - sβ*sε)
	z = R * (ss*sε + sβ*cε)
	return
}

func eqProperMotionToEcl(mα unit.HourAngle, mδ unit.Angle, epoch float64, pos *coord.Ecliptic) (mλ, mβ unit.Angle) {
	ε := nutation.MeanObliquity(base.JulianYearToJDE(epoch))
	sε, cε := ε.Sincos()
	α, δ := coord.EclToEq(pos.Lon, pos.Lat, sε, cε)
	sα, cα := α.Sincos()
	sδ, cδ := δ.Sincos()
	cβ := pos.Lat.Cos()
	mλ = (mδ.Mul(sε*cα) + unit.Angle(mα).Mul(cδ*(cε*cδ+sε*sδ*sα))).Div(cβ * cβ)
	mβ = (mδ.Mul(cε*cδ+sε*sδ*sα) - unit.Angle(mα).Mul(sε*cα*cδ)).Div(cβ)
	return
}

// Nutation returns corrections due to nutation for equatorial coordinates
// of an object.
//
// Results are invalid for objects very near the celestial poles.
func Nutation(α, δ, jd float64) (Δα1, Δδ1 float64) {
	ε := nutation.MeanObliquity(jd)
	sε, cε := math.Sincos(ε)
	Δψ, Δε := nutation.Nutation(jd)
	sα, cα := math.Sincos(α)
	tδ := math.Tan(δ)
	// (23.1) p. 151
	Δα1 = (cε+sε*sα*tδ)*Δψ - cα*tδ*Δε
	Δδ1 = sε*cα*Δψ + sα*Δε
	return
}

func eqProperMotionToEcl(mα, mδ, epoch float64, pos *coord.Ecliptic) (mλ, mβ float64) {
	ε := nutation.MeanObliquity(base.JulianYearToJDE(epoch))
	sε, cε := math.Sincos(ε)
	α, δ := coord.EclToEq(pos.Lon, pos.Lat, sε, cε)
	sα, cα := math.Sincos(α)
	sδ, cδ := math.Sincos(δ)
	cβ := math.Cos(pos.Lat)
	mλ = (mδ*sε*cα + mα*cδ*(cε*cδ+sε*sδ*sα)) / (cβ * cβ)
	mβ = (mδ*(cε*cδ+sε*sδ*sα) - mα*sε*cα*cδ) / cβ
	return
}

// Nutation returns corrections due to nutation for equatorial coordinates
// of an object.
//
// Results are invalid for objects very near the celestial poles.
func Nutation(α unit.RA, δ unit.Angle, jd float64) (Δα1 unit.HourAngle, Δδ1 unit.Angle) {
	ε := nutation.MeanObliquity(jd)
	sε, cε := ε.Sincos()
	Δψ, Δε := nutation.Nutation(jd)
	sα, cα := α.Sincos()
	tδ := δ.Tan()
	// (23.1) p. 151
	Δα1 = unit.HourAngle((cε+sε*sα*tδ)*Δψ.Rad() - cα*tδ*Δε.Rad())
	Δδ1 = Δψ.Mul(sε*cα) + Δε.Mul(sα)
	return
}

// ApparentEquatorialVSOP87 returns the apparent position of the sun as equatorial coordinates.
//
// Result computed by VSOP87, at equator and equinox of date in the FK5 frame,
// and includes effects of nutation and aberration.
//
//	α: right ascension in radians
//	δ: declination in radians
//	R: range in AU
func ApparentEquatorialVSOP87(e *pp.V87Planet, jde float64) (α, δ, R float64) {
	// note: duplicate code from ApparentVSOP87 so we can keep Δε.
	// see also duplicate code in time.E().
	s, β, R := TrueVSOP87(e, jde)
	Δψ, Δε := nutation.Nutation(jde)
	a := aberration(R)
	λ := s + Δψ + a
	ε := nutation.MeanObliquity(jde) + Δε
	sε, cε := math.Sincos(ε)
	α, δ = coord.EclToEq(λ, β, sε, cε)
	return
}

// ApparentEquatorial returns the apparent position of the Sun as equatorial coordinates.
//
//	α: right ascension in radians
//	δ: declination in radians
func ApparentEquatorial(jde float64) (α unit.RA, δ unit.Angle) {
	T := base.J2000Century(jde)
	λ := ApparentLongitude(T)
	ε := nutation.MeanObliquity(jde)
	sλ, cλ := λ.Sincos()
	// (25.8) p. 165
	ε += unit.AngleFromDeg(.00256).Mul(node(T).Cos())
	sε, cε := ε.Sincos()
	α = unit.RAFromRad(math.Atan2(cε*sλ, cλ))
	δ = unit.Angle(math.Asin(sε * sλ))
	return
}

// ESmart computes the "equation of time" for the given JDE.
//
// Result is equation of time as an hour angle.
//
// Result is less accurate that E() but the function has the advantage
// of not requiring the V87Planet object.
func ESmart(jde float64) unit.HourAngle {
	ε := nutation.MeanObliquity(jde)
	t := ε.Mul(.5).Tan()
	y := t * t
	T := base.J2000Century(jde)
	L0 := l0(T * .1)
	e := solar.Eccentricity(T)
	M := solar.MeanAnomaly(T)
	s2L0, c2L0 := L0.Mul(2).Sincos()
	sM := M.Sin()
	// (28.3) p. 185, with double angle identity
	return unit.HourAngle(y*s2L0 - 2*e*sM + 4*e*y*sM*c2L0 -
		y*y*s2L0*c2L0 - 1.25*e*e*M.Mul(2).Sin())
}

// ESmart computes the "equation of time" for the given JDE.
//
// Result is equation of time as an hour angle in radians.
//
// Result is less accurate that E() but the function has the advantage
// of not requiring the V87Planet object.
func ESmart(jde float64) float64 {
	ε := nutation.MeanObliquity(jde)
	t := math.Tan(ε * .5)
	y := t * t
	T := base.J2000Century(jde)
	L0 := l0(T * .1)
	e := solar.Eccentricity(T)
	M := solar.MeanAnomaly(T)
	s2L0, c2L0 := math.Sincos(2 * L0)
	sM := math.Sin(M)
	// (28.3) p. 185
	return y*s2L0 - 2*e*sM + 4*e*y*sM*c2L0 -
		y*y*s2L0*c2L0 - 1.25*e*e*math.Sin(2*M)
}

func ExampleNutation() {
	// Example 22.a, p. 148.
	jd := julian.CalendarGregorianToJD(1987, 4, 10)
	Δψ, Δε := nutation.Nutation(jd)
	ε0 := nutation.MeanObliquity(jd)
	ε := ε0 + Δε
	fmt.Printf("%+.3d\n", sexa.FmtAngle(Δψ))
	fmt.Printf("%+.3d\n", sexa.FmtAngle(Δε))
	fmt.Printf("%.3d\n", sexa.FmtAngle(ε0))
	fmt.Printf("%.3d\n", sexa.FmtAngle(ε))
	// Output:
	// -3″.788
	// +9″.443
	// 23°26′27″.407
	// 23°26′36″.850
}

// 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 in radians.
func E(jde float64, e *pp.V87Planet) float64 {
	τ := 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 := -20.4898 / 3600 * math.Pi / 180 / R
	λ := s + Δψ + a
	ε := nutation.MeanObliquity(jde) + Δε
	sε, cε := math.Sincos(ε)
	α, _ := coord.EclToEq(λ, β, sε, cε)
	// (28.1) p. 183
	E := L0 - .0057183*math.Pi/180 - α + Δψ*cε
	return base.PMod(E+math.Pi, 2*math.Pi) - math.Pi
}

// 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)
}

// Aberration returns corrections due to aberration for equatorial
// coordinates of an object.
func Aberration(α, δ, jd float64) (Δα2, Δδ2 float64) {
	ε := nutation.MeanObliquity(jd)
	T := base.J2000Century(jd)
	s, _ := solar.True(T)
	e := solar.Eccentricity(T)
	π := perihelion(T)
	sα, cα := math.Sincos(α)
	sδ, cδ := math.Sincos(δ)
	ss, cs := math.Sincos(s)
	sπ, cπ := math.Sincos(π)
	cε := math.Cos(ε)
	tε := math.Tan(ε)
	q1 := cα * cε
	// (23.3) p. 152
	Δα2 = κ * (e*(q1*cπ+sα*sπ) - (q1*cs + sα*ss)) / cδ
	q2 := cε * (tε*cδ - sα*sδ)
	q3 := cα * sδ
	Δδ2 = κ * (e*(cπ*q2+sπ*q3) - (cs*q2 + ss*q3))
	return
}

// Aberration returns corrections due to aberration for equatorial
// coordinates of an object.
func Aberration(α unit.RA, δ unit.Angle, jd float64) (Δα2 unit.HourAngle, Δδ2 unit.Angle) {
	ε := nutation.MeanObliquity(jd)
	T := base.J2000Century(jd)
	s, _ := solar.True(T)
	e := solar.Eccentricity(T)
	π := perihelion(T)
	sα, cα := α.Sincos()
	sδ, cδ := δ.Sincos()
	ss, cs := s.Sincos()
	sπ, cπ := π.Sincos()
	cε := ε.Cos()
	tε := ε.Tan()
	q1 := cα * cε
	// (23.3) p. 152
	Δα2 = unit.HourAngle(κ.Rad() * (e*(q1*cπ+sα*sπ) - (q1*cs + sα*ss)) / cδ)
	q2 := cε * (tε*cδ - sα*sδ)
	q3 := cα * sδ
	Δδ2 = κ.Mul(e*(cπ*q2+sπ*q3) - (cs*q2 + ss*q3))
	return
}

// 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β")
	}
}