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:   11.94524 * math.Pi / 180,
		Node:  334.75006 * math.Pi / 180,
		ArgP:  186.23352 * math.Pi / 180,
	}
	j := julian.CalendarGregorianToJD(1990, 10, 6)
	α, δ, ψ := k.Position(j, earth)
	fmt.Printf("α = %.1d\n", sexa.NewFmtRA(α))
	fmt.Printf("δ = %.0d\n", sexa.NewFmtAngle(δ))
	fmt.Printf("ψ = %.2f\n", ψ*180/math.Pi)
	// Output:
	// α = 10ʰ34ᵐ14ˢ.2
	// δ = 19°9′31″
	// ψ = 40.51
}

func ExampleTime() {
	// Example 19.a, p. 121.
	r1 := unit.AngleFromDeg(113.56833)
	d1 := unit.AngleFromDeg(31.89756)
	r2 := unit.AngleFromDeg(116.25042)
	d2 := unit.AngleFromDeg(28.03681)
	r3 := make([]unit.Angle, 5)
	for i, ri := range []float64{
		118.98067, 119.59396, 120.20413, 120.81108, 121.41475} {
		r3[i] = unit.AngleFromDeg(ri)
	}
	d3 := make([]unit.Angle, 5)
	for i, di := range []float64{
		21.68417, 21.58983, 21.49394, 21.39653, 21.29761} {
		d3[i] = unit.AngleFromDeg(di)
	}
	// use JD as time to handle month boundary
	jd, err := line.Time(r1, d1, r2, d2, r3, d3,
		julian.CalendarGregorianToJD(1994, 9, 29),
		julian.CalendarGregorianToJD(1994, 10, 3))
	if err != nil {
		fmt.Println(err)
		return
	}
	y, m, d := julian.JDToCalendar(jd)
	dInt, dFrac := math.Modf(d)
	fmt.Printf("%d %s %.4f\n", y, time.Month(m), d)
	fmt.Printf("%d %s %d, at %h TD(UT)\n", y, time.Month(m), int(dInt),
		sexa.FmtTime(unit.TimeFromDay(dFrac)))
	// Output:
	// 1994 October 1.2233
	// 1994 October 1, at 5ʰ TD(UT)
}

func ExampleTime() {
	// Example 19.a, p. 121.

	// convert degree data to radians
	r1 := 113.56833 * math.Pi / 180
	d1 := 31.89756 * math.Pi / 180
	r2 := 116.25042 * math.Pi / 180
	d2 := 28.03681 * math.Pi / 180
	r3 := make([]float64, 5)
	for i, ri := range []float64{
		118.98067, 119.59396, 120.20413, 120.81108, 121.41475} {
		r3[i] = ri * math.Pi / 180
	}
	d3 := make([]float64, 5)
	for i, di := range []float64{
		21.68417, 21.58983, 21.49394, 21.39653, 21.29761} {
		d3[i] = di * math.Pi / 180
	}
	// use JD as time to handle month boundary
	jd, err := line.Time(r1, d1, r2, d2, r3, d3,
		julian.CalendarGregorianToJD(1994, 9, 29),
		julian.CalendarGregorianToJD(1994, 10, 3))
	if err != nil {
		fmt.Println(err)
		return
	}
	y, m, d := julian.JDToCalendar(jd)
	dInt, dFrac := math.Modf(d)
	fmt.Printf("%d %s %.4f\n", y, time.Month(m), d)
	fmt.Printf("%d %s %d, at %h TD(UT)\n", y, time.Month(m), int(dInt),
		sexa.NewFmtTime(dFrac*24*3600))
	// Output:
	// 1994 October 1.2233
	// 1994 October 1, at 5ʰ TD(UT)
}

func ExampleCalendarGregorianToJD_halley() {
	// Example 7.c, p. 64.
	jd1 := julian.CalendarGregorianToJD(1910, 4, 20)
	jd2 := julian.CalendarGregorianToJD(1986, 2, 9)
	fmt.Printf("%.0f days\n", jd2-jd1)
	// Output:
	// 27689 days
}

func ExampleElements_AnomalyDistance() {
	// Example 34.a, p. 243
	e := &parabolic.Elements{
		TimeP: julian.CalendarGregorianToJD(1998, 4, 14.4358),
		PDis:  1.487469,
	}
	j := julian.CalendarGregorianToJD(1998, 8, 5)
	ν, r := e.AnomalyDistance(j)
	fmt.Printf("%.5f deg\n", ν*180/math.Pi)
	fmt.Printf("%.6f AU\n", r)
	// Output:
	// 66.78862 deg
	// 2.133911 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″
}

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
}

func ExampleInterp10A() {
	// Example 10.a, p. 78.
	dt := deltat.Interp10A(julian.CalendarGregorianToJD(1977, 2, 18))
	fmt.Printf("%+.1f seconds\n", dt)
	// Output:
	// +47.6 seconds
}

func ExampleV87Planet_Position() {
	// Example 32.a, p. 219
	jd := julian.CalendarGregorianToJD(1992, 12, 20)
	p, err := pp.LoadPlanet(pp.Venus, "")
	if err != nil {
		fmt.Println(err)
		return
	}
	l, b, r := p.Position(jd)
	fmt.Printf("L = %s\n",
		base.DecSymAdd(fmt.Sprintf("%+.5f", l*180/math.Pi), '°'))
	fmt.Printf("B = %s\n",
		base.DecSymAdd(fmt.Sprintf("%+.5f", b*180/math.Pi), '°'))
	fmt.Printf("R = %.6f AU\n", r)
	// Meeus results:
	// L = +26°.11428
	// B = -2°.62070
	// R = 0.724603 AU
	// Answers below seem close enough.

	// Output:
	// L = +26°.11412
	// B = -2°.62060
	// R = 0.724602 AU
}

func ExampleCalendarGregorianToJD_sputnik() {
	// Example 7.a, p. 61.
	jd := julian.CalendarGregorianToJD(1957, 10, 4.81)
	fmt.Printf("%.2f\n", jd)
	// Output:
	// 2436116.31
}

func ExampleApproxTimes() {
	// Example 15.a, p. 103.
	jd := julian.CalendarGregorianToJD(1988, 3, 20)
	p := globe.Coord{
		Lon: base.NewAngle(false, 71, 5, 0).Rad(),
		Lat: base.NewAngle(false, 42, 20, 0).Rad(),
	}
	// Meeus gives us the value of 11h 50m 58.1s but we have a package
	// function for this:
	Th0 := sidereal.Apparent0UT(jd)
	α := base.NewRA(2, 46, 55.51).Rad()
	δ := base.NewAngle(false, 18, 26, 27.3).Rad()
	h0 := rise.Stdh0Stellar
	rise, transit, set, err := rise.ApproxTimes(p, h0, Th0, α, δ)
	if err != nil {
		fmt.Println(err)
		return
	}
	// Units for approximate values given near top of p. 104 are circles.
	fmt.Printf("rising:  %+.5f\n", rise/86400)
	fmt.Printf("transit: %+.5f\n", transit/86400)
	fmt.Printf("seting:  %+.5f\n", set/86400)
	// Output:
	// rising:  +0.51816
	// transit: +0.81965
	// seting:  +0.12113
}

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 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 ExampleParallax() {
	// Example 47.a, p. 342.
	_, _, Δ := moon.Position(julian.CalendarGregorianToJD(1992, 4, 12))
	π := moon.Parallax(Δ)
	fmt.Printf("π = %.6f\n", π*180/math.Pi)
	// Output:
	// π = 0.991990
}

func TestNode(t *testing.T) {
	j := julian.CalendarGregorianToJD(2065, 6, 24)
	var e pe.Elements
	pe.Mean(pe.Mercury, j, &e)
	if Ω := pe.Node(pe.Mercury, j); Ω != e.Node {
		t.Fatal(Ω, "!=", e.Node)
	}
}

func TestInc(t *testing.T) {
	j := julian.CalendarGregorianToJD(2065, 6, 24)
	var e pe.Elements
	pe.Mean(pe.Mercury, j, &e)
	if i := pe.Inc(pe.Mercury, j); i != e.Inc {
		t.Fatal(i, "!=", e.Inc)
	}
}

func ExampleESmart() {
	// Example 28.b, p. 185
	eq := eqtime.ESmart(julian.CalendarGregorianToJD(1992, 10, 13))
	fmt.Printf("+%.7f rad\n", eq)
	fmt.Printf("%+.1d", sexa.FmtHourAngle(eq))
	// Output:
	// +0.0598256 rad
	// +13ᵐ42ˢ.7
}

func ExamplePhaseAngle3() {
	i := moonillum.PhaseAngle3(julian.CalendarGregorianToJD(1992, 4, 12))
	k := base.Illuminated(i)
	fmt.Printf("i = %.2f\n", i*180/math.Pi)
	fmt.Printf("k = %.4f\n", k)
	// Output:
	// i = 68.88
	// k = 0.6801
}

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

// Planet computes UT rise, transit and set times for a planet on a day of
// interest.
//
//  yr, mon, day are the Gregorian date.
//  pos is geographic coordinates of observer.
//  e must be a V87Planet object for Earth
//  pl must be a V87Planet object for another planet.
//
// Obtain V87Planet objects with the planetposition package.
//
// Result units are seconds of day and are in the range [0,86400).
func Planet(yr, mon, day int, pos globe.Coord, e, pl *pp.V87Planet) (tRise, tTransit, tSet unit.Time, err error) {
	jd := julian.CalendarGregorianToJD(yr, mon, float64(day))
	α := make([]unit.RA, 3)
	δ := make([]unit.Angle, 3)
	α[0], δ[0] = elliptic.Position(pl, e, jd-1)
	α[1], δ[1] = elliptic.Position(pl, e, jd)
	α[2], δ[2] = elliptic.Position(pl, e, jd+1)
	return Times(pos, deltat.Interp10A(jd), Stdh0Stellar,
		sidereal.Apparent0UT(jd), α, δ)
}