func ExampleGeneral_b() {
	// Example 58.b, p. 404.
	const p = math.Pi / 180
	ls, c, _, ψ := sundial.General(
		unit.AngleFromDeg(-35),
		unit.AngleFromDeg(160),
		1,
		unit.AngleFromDeg(90))
	for _, l := range ls {
		if l.Hour == 12 {
			fmt.Printf("%d:  x = %+.4f  y = %+.4f\n",
				l.Hour, l.Points[5].X, l.Points[5].Y)
		}
		if l.Hour == 15 {
			fmt.Printf("%d:  x = %+.4f  y = %+.4f\n",
				l.Hour, l.Points[3].X, l.Points[3].Y)
		}
	}
	fmt.Printf("x0 = %+.4f\n", c.X)
	fmt.Printf("y0 = %+.4f\n", c.Y)
	fmt.Printf("ψ = %.4f\n", ψ.Deg())
	// Output:
	// 12:  x = +0.3640  y = -0.7410
	// 15:  x = -0.8439  y = -0.9298
	// x0 = +0.3640
	// y0 = +0.7451
	// ψ = 50.3315
}

func ExampleApparentEccentricity() {
	// Example 57.b, p. 400
	fmt.Printf("%.3f\n", binary.ApparentEccentricity(.2763,
		unit.AngleFromDeg(59.025), unit.AngleFromDeg(219.907)))
	// Output:
	// 0.860
}

func ExampleGeneral_a() {
	// Example 58.a, p. 404.
	ls, c, _, ψ := sundial.General(
		unit.AngleFromDeg(40),
		unit.AngleFromDeg(70),
		1,
		unit.AngleFromDeg(50))
	fmt.Printf("Hours:  %d", ls[0].Hour)
	for _, l := range ls[1:] {
		fmt.Printf(", %d", l.Hour)
	}
	fmt.Println()
	for _, l := range ls {
		if l.Hour == 11 {
			fmt.Printf("%d:  x = %.4f  y = %.4f\n",
				l.Hour, l.Points[2].X, l.Points[2].Y)
		}
		if l.Hour == 14 {
			fmt.Printf("%d:  x = %.4f  y = %.4f\n",
				l.Hour, l.Points[6].X, l.Points[6].Y)
		}
	}
	fmt.Printf("x0 = %+.4f\n", c.X)
	fmt.Printf("y0 = %+.4f\n", c.Y)
	fmt.Printf("ψ = %.4f\n", ψ.Deg())
	// Output:
	// Hours:  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
	// 11:  x = -2.0007  y = -1.1069
	// 14:  x = -0.0390  y = -0.3615
	// x0 = +3.3880
	// y0 = -3.1102
	// ψ = 12.2672
}

func ExampleElements_Position() {
	// Example 33.b, p. 232.
	earth, err := pp.LoadPlanet(pp.Earth)
	if err != nil {
		fmt.Println(err)
		return
	}
	k := &elliptic.Elements{
		TimeP: julian.CalendarGregorianToJD(1990, 10, 28.54502),
		Axis:  2.2091404,
		Ecc:   .8502196,
		Inc:   unit.AngleFromDeg(11.94524),
		Node:  unit.AngleFromDeg(334.75006),
		ArgP:  unit.AngleFromDeg(186.23352),
	}
	j := julian.CalendarGregorianToJD(1990, 10, 6)
	α, δ, ψ := k.Position(j, earth)
	fmt.Printf("α = %.1d\n", sexa.FmtRA(α))
	fmt.Printf("δ = %.0d\n", sexa.FmtAngle(δ))
	fmt.Printf("ψ = %.2f\n", ψ.Deg())
	// Output:
	// α = 10ʰ34ᵐ14ˢ.2
	// δ = 19°9′31″
	// ψ = 40.51
}

func ExampleSaturn() {
	// Example 41.d, p. 285
	v := illum.Saturn(9.867882, 10.464606,
		unit.AngleFromDeg(16.442), unit.AngleFromDeg(4.198))
	fmt.Printf("%+.1f\n", v)
	// Output:
	// +0.9
}

func ExamplePhaseAngle3() {
	// Example 41.a, p. 284
	i := illum.PhaseAngle3(
		unit.AngleFromDeg(26.10588),
		unit.AngleFromDeg(-2.62102),
		.621794, -.664905, -.033138, .910947)
	fmt.Printf("%.5f\n", i.Cos())
	// Output:
	// 0.29312
}

// True returns true geometric longitude and anomaly of the sun referenced to the mean equinox of date.
//
// Argument T is the number of Julian centuries since J2000.
// See base.J2000Century.
//
// Results:
//	s = true geometric longitude, ☉
//	ν = true anomaly
func True(T float64) (s, ν unit.Angle) {
	// (25.2) p. 163
	L0 := unit.AngleFromDeg(base.Horner(T, 280.46646, 36000.76983, 0.0003032))
	M := MeanAnomaly(T)
	C := unit.AngleFromDeg(base.Horner(T, 1.914602, -0.004817, -.000014)*
		M.Sin() +
		(0.019993-.000101*T)*M.Mul(2).Sin() +
		0.000289*M.Mul(3).Sin())
	return (L0 + C).Mod1(), (M + C).Mod1()
}

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

func ExampleLimb() {
	// Example 48.a, p. 347.
	χ := base.Limb(
		unit.RAFromDeg(134.6885),
		unit.AngleFromDeg(13.7684),
		unit.RAFromDeg(20.6579),
		unit.AngleFromDeg(8.6964))
	fmt.Printf("χ = %.1f\n", χ.Deg())
	// Output:
	// χ = 285.0
}

func ExamplePhaseAngleEq2() {
	i := moonillum.PhaseAngleEq2(
		unit.RAFromDeg(134.6885),
		unit.AngleFromDeg(13.7684),
		unit.RAFromDeg(20.6579),
		unit.AngleFromDeg(8.6964))
	k := base.Illuminated(i)
	fmt.Printf("k = %.4f\n", k)
	// Output:
	// k = 0.6775
}

func ExamplePhaseAngleEq() {
	i := moonillum.PhaseAngleEq(
		unit.RAFromDeg(134.6885),
		unit.AngleFromDeg(13.7684),
		368410,
		unit.RAFromDeg(20.6579),
		unit.AngleFromDeg(8.6964),
		149971520)
	fmt.Printf("i = %.4f\n", i.Deg())
	// Output:
	// i = 69.0756
}

func ExamplePhaseAngle2() {
	// Example 41.a, p. 284
	i := illum.PhaseAngle2(
		unit.AngleFromDeg(26.10588),
		unit.AngleFromDeg(-2.62102),
		.724604,
		unit.AngleFromDeg(88.35704),
		.983824, .910947)
	fmt.Printf("%.5f\n", i.Cos())
	// Output:
	// 0.29312
}

func ExampleEclipticAtHorizon() {
	ε := unit.AngleFromDeg(23.44)
	φ := unit.AngleFromDeg(51)
	θ := unit.TimeFromHour(5)
	λ1, λ2, I := parallactic.EclipticAtHorizon(ε, φ, θ)
	fmt.Println(sexa.FmtAngle(λ1))
	fmt.Println(sexa.FmtAngle(λ2))
	fmt.Println(sexa.FmtAngle(I))
	// Output:
	// 169°21′30″
	// 349°21′30″
	// 61°53′14″
}

func ExampleSunAltitude() {
	j := julian.CalendarGregorianToJD(1992, 4, 12)
	earth, err := pp.LoadPlanet(pp.Earth)
	if err != nil {
		fmt.Println(err)
		return
	}
	_, _, _, l0, b0 := moon.Physical(j, earth)
	h := moon.SunAltitude(
		unit.AngleFromDeg(-20), unit.AngleFromDeg(9.7), l0, b0)
	fmt.Printf("%+.3f\n", h.Deg())
	// Output:
	// +2.318
}

func ExampleSunrise() {
	earth, err := pp.LoadPlanet(pp.Earth)
	if err != nil {
		fmt.Println(err)
		return
	}
	j0 := julian.CalendarGregorianToJD(1992, 4, 15)
	j := moon.Sunrise(
		unit.AngleFromDeg(-20), unit.AngleFromDeg(9.7), j0, earth)
	y, m, d := julian.JDToCalendar(j)
	fmt.Printf("%d %s %.4f TD\n", y, time.Month(m), d)
	// Output:
	// 1992 April 11.8069 TD
}

func TestDiurnalPathAtHorizon(t *testing.T) {
	φ := unit.AngleFromDeg(40)
	ε := unit.AngleFromDeg(23.44)
	J := parallactic.DiurnalPathAtHorizon(0, φ)
	Jexp := math.Pi/2 - φ
	if math.Abs((J-Jexp).Rad()/Jexp.Rad()) > 1e-15 {
		t.Fatal("0 dec:", sexa.FmtAngle(J))
	}
	J = parallactic.DiurnalPathAtHorizon(ε, φ)
	Jexp = unit.NewAngle(' ', 45, 31, 0)
	if math.Abs((J-Jexp).Rad()/Jexp.Rad()) > 1e-3 {
		t.Fatal("solstice:", sexa.FmtAngle(J))
	}
}

func ExampleGeneral_c() {
	// Example 58.c, p. 405.
	ls, _, _, _ := sundial.General(
		unit.AngleFromDeg(40),
		unit.AngleFromDeg(160),
		1,
		unit.AngleFromDeg(75))
	fmt.Printf("Hours:  %d", ls[0].Hour)
	for _, l := range ls[1:] {
		fmt.Printf(", %d", l.Hour)
	}
	fmt.Println()
	// Output:
	// Hours:  5, 6, 13, 14, 15, 16, 17, 18, 19
}

func ExampleEclipticPosition() {
	// Example 21.c, p. 137.
	eclFrom := &coord.Ecliptic{
		Lat: unit.AngleFromDeg(1.76549),
		Lon: unit.AngleFromDeg(149.48194),
	}
	eclTo := &coord.Ecliptic{}
	epochFrom := 2000.0
	epochTo := base.JDEToJulianYear(julian.CalendarJulianToJD(-214, 6, 30))
	precess.EclipticPosition(eclFrom, eclTo, epochFrom, epochTo, 0, 0)
	fmt.Printf("%.3f\n", eclTo.Lon.Deg())
	fmt.Printf("%+.3f\n", eclTo.Lat.Deg())
	// Output:
	// 118.704
	// +1.615
}

// PhaseAngle3 computes the phase angle of the Moon given a julian day.
//
// Less accurate than PhaseAngle functions taking coordinates.
func PhaseAngle3(jde float64) unit.Angle {
	T := base.J2000Century(jde)
	D := unit.AngleFromDeg(base.Horner(T, 297.8501921,
		445267.1114034, -.0018819, 1/545868, -1/113065000)).Mod1().Rad()
	M := unit.AngleFromDeg(base.Horner(T,
		357.5291092, 35999.0502909, -.0001535, 1/24490000)).Mod1().Rad()
	Mʹ := unit.AngleFromDeg(base.Horner(T, 134.9633964,
		477198.8675055, .0087414, 1/69699, -1/14712000)).Mod1().Rad()
	return math.Pi - unit.Angle(D) + unit.AngleFromDeg(
		-6.289*math.Sin(Mʹ)+
			2.1*math.Sin(M)+
			-1.274*math.Sin(2*D-Mʹ)+
			-.658*math.Sin(2*D)+
			-.214*math.Sin(2*Mʹ)+
			-.11*math.Sin(D))
}

func ExampleVenus() {
	// Example 41.c, p. 285
	v := illum.Venus(.724604, .910947, unit.AngleFromDeg(72.96))
	fmt.Printf("%.1f\n", v)
	// Output:
	// -3.8
}