func ExampleEccentricity() {
	// Example 25.a, p. 165.
	T := base.J2000Century(julian.CalendarGregorianToJD(1992, 10, 13))
	fmt.Printf("%.9f\n", solar.Eccentricity(T))
	// Output:
	// 0.016711668
}

func ExampleMeanAnomaly() {
	// Example 25.a, p. 165.
	T := base.J2000Century(julian.CalendarGregorianToJD(1992, 10, 13))
	fmt.Printf("%.5f\n", solar.MeanAnomaly(T)*180/math.Pi)
	// Output:
	// -2241.00603
}

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

func ExampleRadius() {
	// Example 25.a, p. 165.
	T := base.J2000Century(julian.CalendarGregorianToJD(1992, 10, 13))
	fmt.Printf("%.5f AU\n", solar.Radius(T))
	// Output:
	// 0.99766 AU
}

func ExampleApparentLongitude() {
	// Example 25.a, p. 165.
	T := base.J2000Century(julian.CalendarGregorianToJD(1992, 10, 13))
	fmt.Println("λ:", sexa.NewFmtAngle(solar.ApparentLongitude(T)))
	// Output:
	// λ: 199°54′32″
}

// AberrationRonVondrak uses the Ron-Vondrák expression to compute corrections
// due to aberration for equatorial coordinates of an object.
func AberrationRonVondrak(α, δ, jd float64) (Δα, Δδ float64) {
	T := base.J2000Century(jd)
	r := &rv{
		T:  T,
		L2: 3.1761467 + 1021.3285546*T,
		L3: 1.7534703 + 628.3075849*T,
		L4: 6.2034809 + 334.0612431*T,
		L5: 0.5995465 + 52.9690965*T,
		L6: 0.8740168 + 21.3299095*T,
		L7: 5.4812939 + 7.4781599*T,
		L8: 5.3118863 + 3.8133036*T,
		Lp: 3.8103444 + 8399.6847337*T,
		D:  5.1984667 + 7771.3771486*T,
		Mp: 2.3555559 + 8328.6914289*T,
		F:  1.6279052 + 8433.4661601*T,
	}
	var Xp, Yp, Zp float64
	// sum smaller terms first
	for i := 35; i >= 0; i-- {
		x, y, z := rvTerm[i](r)
		Xp += x
		Yp += y
		Zp += z
	}
	sα, cα := math.Sincos(α)
	sδ, cδ := math.Sincos(δ)
	// (23.4) p. 156
	return (Yp*cα - Xp*sα) / (c * cδ), -((Xp*cα+Yp*sα)*sδ - Zp*cδ) / c
}

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

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

// 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.
//
// Result unit is radians.
func MeanObliquity(jde float64) float64 {
	// (22.2) p. 147
	return base.Horner(base.J2000Century(jde),
		base.NewAngle(false, 23, 26, 21.448).Rad(),
		-46.815/3600*(math.Pi/180),
		-0.00059/3600*(math.Pi/180),
		0.001813/3600*(math.Pi/180))
}

// Mean computes some intermediate values for a mean planetary configuration
// given a year and a row of coefficients from Table 36.A, p. 250.
func mean(y float64, a *ca) (J, M, T float64) {
	// (36.1) p. 250
	k := math.Floor((365.2425*y+1721060-a.A)/a.B + .5)
	J = a.A + k*a.B
	M = unit.PMod(a.M0+k*a.M1, 360) * math.Pi / 180
	T = base.J2000Century(J)
	return
}

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

// TrueNode returns longitude of the true ascending node.
//
// That is, the node of the instantaneous lunar orbit.
//
// Result in radians.
func TrueNode(jde float64) float64 {
	D, M, Mʹ, F := dmf(base.J2000Century(jde))
	return Node(jde) +
		-1.4979*p*math.Sin(2*(D-F)) +
		-.15*p*math.Sin(M) +
		-.1226*p*math.Sin(2*D) +
		.1176*p*math.Sin(2*F) +
		-.0801*p*math.Sin(2*(Mʹ-F))
}

// TrueNode returns longitude of the true ascending node.
//
// That is, the node of the instantaneous lunar orbit.
func TrueNode(jde float64) unit.Angle {
	D, M, Mʹ, F := dmf(base.J2000Century(jde))
	return Node(jde) + unit.AngleFromDeg(
		-1.4979*math.Sin(2*(D-F))+
			-.15*math.Sin(M)+
			-.1226*math.Sin(2*D)+
			.1176*math.Sin(2*F)+
			-.0801*math.Sin(2*(Mʹ-F)))
}

// cl splits the work into two closures.
func cl(jde float64, earth, saturn *pp.V87Planet) (f1 func() (ΔU, B float64),
	f2 func() (Bʹ, P, aEdge, bEdge float64)) {
	const p = math.Pi / 180
	var i, Ω float64
	var l0, b0, R float64
	Δ := 9.
	var λ, β float64
	var si, ci, sβ, cβ, sB float64
	var sbʹ, cbʹ, slʹΩ, clʹΩ float64
	f1 = func() (ΔU, B float64) {
		// (45.1), p. 318
		T := base.J2000Century(jde)
		i = base.Horner(T, 28.075216*p, -.012998*p, .000004*p)
		Ω = base.Horner(T, 169.50847*p, 1.394681*p, .000412*p)
		// Step 2.
		l0, b0, R = earth.Position(jde)
		l0, b0 = pp.ToFK5(l0, b0, jde)
		sl0, cl0 := math.Sincos(l0)
		sb0 := math.Sin(b0)
		// Steps 3, 4.
		var l, b, r, x, y, z float64
		f := func() {
			τ := base.LightTime(Δ)
			l, b, r = saturn.Position(jde - τ)
			l, b = pp.ToFK5(l, b, jde)
			sl, cl := math.Sincos(l)
			sb, cb := math.Sincos(b)
			x = r*cb*cl - R*cl0
			y = r*cb*sl - R*sl0
			z = r*sb - R*sb0
			Δ = math.Sqrt(x*x + y*y + z*z)
		}
		f()
		f()
		// Step 5.
		λ = math.Atan2(y, x)
		β = math.Atan(z / math.Hypot(x, y))
		// First part of step 6.
		si, ci = math.Sincos(i)
		sβ, cβ = math.Sincos(β)
		sB = si*cβ*math.Sin(λ-Ω) - ci*sβ
		B = math.Asin(sB) // return value
		// Step 7.
		N := 113.6655*p + .8771*p*T
		lʹ := l - .01759*p/r
		bʹ := b - .000764*p*math.Cos(l-N)/r
		// Setup for steps 8, 9.
		sbʹ, cbʹ = math.Sincos(bʹ)
		slʹΩ, clʹΩ = math.Sincos(lʹ - Ω)
		// Step 9.
		sλΩ, cλΩ := math.Sincos(λ - Ω)
		U1 := math.Atan2(si*sbʹ+ci*cbʹ*slʹΩ, cbʹ*clʹΩ)
		U2 := math.Atan2(si*sβ+ci*cβ*sλΩ, cβ*cλΩ)
		ΔU = math.Abs(U1 - U2) // return value
		return
	}

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

// Mean returns mean orbital elements for a planet
//
// Argument p must be a planet const as defined above, argument e is
// a result parameter.  A valid non-nil pointer to an Elements struct
// must be passed in.
//
// Results are referenced to mean dynamical ecliptic and equinox of date.
//
// Semimajor axis is in AU, angular elements are in radians.
func Mean(p int, jde float64, e *Elements) {
	T := base.J2000Century(jde)
	c := &cMean[p]
	e.Lon = unit.AngleFromDeg(base.Horner(T, c.L...)).Mod1()
	e.Axis = base.Horner(T, c.a...)
	e.Ecc = base.Horner(T, c.e...)
	e.Inc = unit.AngleFromDeg(base.Horner(T, c.i...))
	e.Node = unit.AngleFromDeg(base.Horner(T, c.Ω...))
	e.Peri = unit.AngleFromDeg(base.Horner(T, c.ϖ...))
}

// Mean returns mean orbital elements for a planet
//
// Argument p must be a planet const as defined above, argument e is
// a result parameter.  A valid non-nil pointer to an Elements struct
// must be passed in.
//
// Results are referenced to mean dynamical ecliptic and equinox of date.
//
// Semimajor axis is in AU, angular elements are in radians.
func Mean(p int, jde float64, e *Elements) {
	T := base.J2000Century(jde)
	c := &cMean[p]
	e.Lon = base.PMod(base.Horner(T, c.L...)*math.Pi/180, 2*math.Pi)
	e.Axis = base.Horner(T, c.a...)
	e.Ecc = base.Horner(T, c.e...)
	e.Inc = base.Horner(T, c.i...) * math.Pi / 180
	e.Node = base.Horner(T, c.Ω...) * math.Pi / 180
	e.Peri = base.Horner(T, c.ϖ...) * math.Pi / 180
}

// FractionVenus computes an approximation of the illumanted fraction of Venus.
func FractionVenus(jde float64) float64 {
	T := base.J2000Century(jde)
	V := 261.51*p + 22518.443*p*T
	M := 177.53*p + 35999.05*p*T
	N := 50.42*p + 58517.811*p*T
	W := V + 1.91*p*math.Sin(M) + .78*p*math.Sin(N)
	Δ := math.Sqrt(1.52321 + 1.44666*math.Cos(W))
	s := .72333 + Δ
	return (s*s - 1) / 2.89332 / Δ
}

// 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 in radians
//	β: ecliptic latitude in radians
//	R: range in AU
func TrueVSOP87(e *pp.V87Planet, jde float64) (s, β, R float64) {
	l, b, r := e.Position(jde)
	s = l + math.Pi
	// FK5 correction.
	λp := base.Horner(base.J2000Century(jde),
		s, -1.397*math.Pi/180, -.00031*math.Pi/180)
	sλp, cλp := math.Sincos(λp)
	Δβ := .03916 / 3600 * math.Pi / 180 * (cλp - sλp)
	// (25.9) p. 166
	return base.PMod(s-.09033/3600*math.Pi/180, 2*math.Pi), Δβ - b, r
}

// ToFK5 converts ecliptic longitude and latitude from dynamical frame to FK5.
func ToFK5(L, B, jde float64) (L5, B5 float64) {
	// formula 32.3, p. 219.
	T := base.J2000Century(jde)
	Lp := L - 1.397*math.Pi/180*T - .00031*math.Pi/180*T*T
	sLp, cLp := math.Sincos(Lp)
	// (32.3) p. 219
	L5 = L + -.09033/3600*math.Pi/180 +
		.03916/3600*math.Pi/180*(cLp+sLp)*math.Tan(B)
	B5 = B + .03916/3600*math.Pi/180*(cLp-sLp)
	return
}