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

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

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

func xyz(e *pp.V87Planet, jde float64) (x, y, z float64) {
	l, b, r := e.Position2000(jde)
	s := l + math.Pi
	β := -b
	ss, cs := math.Sincos(s)
	sβ, cβ := math.Sincos(β)
	// (26.2) p. 172
	x = r * cβ * cs
	y = r * cβ * ss
	z = r * sβ
	return
}

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

func ap2a(j1, d float64, a bool, v *pp.V87Planet) (jde, r float64) {
	// Meeus doesn't give a complete algorithm for finding accurate answers.
	// The algorithm here starts with the approximate result algorithm
	// then crawls along the curve at resultion d until three points
	// are found containing the desired extremum.  It's slow if the starting
	// point is far away (the case of Neptune) or if high accuracy is
	// demanded.  1 day accuracy is generally quick.
	j0 := j1 - d
	j2 := j1 + d
	rr := make([]float64, 3)
	_, _, rr[0] = v.Position2000(j0)
	_, _, rr[1] = v.Position2000(j1)
	_, _, rr[2] = v.Position2000(j2)
	for {
		if a {
			if rr[1] > rr[0] && rr[1] > rr[2] {
				break
			}
		} else {
			if rr[1] < rr[0] && rr[1] < rr[2] {
				break
			}
		}
		if rr[0] < rr[2] == a {
			j0 = j1
			j1 = j2
			j2 += d
			rr[0] = rr[1]
			rr[1] = rr[2]
			_, _, rr[2] = v.Position2000(j2)
		} else {
			j2 = j1
			j1 = j0
			j0 -= d
			rr[2] = rr[1]
			rr[1] = rr[0]
			_, _, rr[0] = v.Position2000(j0)
		}
	}
	l, err := interp.NewLen3(j0, j2, rr)
	if err != nil {
		panic(err) // unexpected.
	}
	jde, r, err = l.Extremum()
	if err != nil {
		panic(err) // unexpected.
	}
	return
}

// Positions returns positions of the eight major moons of Saturn.
//
// Results returned in argument pos, which must not be nil.
//
// Result units are Saturn radii.
func Positions(jde float64, earth, saturn *pp.V87Planet, pos *[8]XY) {
	s, β, R := solar.TrueVSOP87(earth, jde)
	ss, cs := s.Sincos()
	sβ := β.Sin()
	Δ := 9.
	var x, y, z float64
	var JDE float64
	f := func() {
		τ := base.LightTime(Δ)
		JDE = jde - τ
		l, b, r := saturn.Position(JDE)
		l, b = pp.ToFK5(l, b, JDE)
		sl, cl := l.Sincos()
		sb, cb := b.Sincos()
		x = r*cb*cl + R*cs
		y = r*cb*sl + R*ss
		z = r*sb + R*sβ
		Δ = math.Sqrt(x*x + y*y + z*z)
	}
	f()
	f()
	λ0 := unit.Angle(math.Atan2(y, x))
	β0 := unit.Angle(math.Atan(z / math.Hypot(x, y)))
	ecl := &coord.Ecliptic{λ0, β0}
	precess.EclipticPosition(ecl, ecl,
		base.JDEToJulianYear(jde), base.JDEToJulianYear(base.B1950), 0, 0)
	λ0, β0 = ecl.Lon, ecl.Lat
	q := newQs(JDE)
	s4 := [9]r4{{}, // 0 unused
		q.mimas(),
		q.enceladus(),
		q.tethys(),
		q.dione(),
		q.rhea(),
		q.titan(),
		q.hyperion(),
		q.iapetus(),
	}
	var X, Y, Z [9]float64
	for j := 1; j <= 8; j++ {
		u := s4[j].λ - s4[j].Ω
		w := s4[j].Ω - 168.8112*d
		su, cu := math.Sincos(u)
		sw, cw := math.Sincos(w)
		sγ, cγ := math.Sincos(s4[j].γ)
		r := s4[j].r
		X[j] = r * (cu*cw - su*cγ*sw)
		Y[j] = r * (su*cw*cγ + cu*sw)
		Z[j] = r * su * sγ
	}
	Z[0] = 1
	sλ0, cλ0 := λ0.Sincos()
	sβ0, cβ0 := β0.Sincos()
	var A, B, C [9]float64
	for j := range X {
		a := X[j]
		b := q.c1*Y[j] - q.s1*Z[j]
		c := q.s1*Y[j] + q.c1*Z[j]
		a, b =
			q.c2*a-q.s2*b,
			q.s2*a+q.c2*b
		A[j], b =
			a*sλ0-b*cλ0,
			a*cλ0+b*sλ0
		B[j], C[j] =
			b*cβ0+c*sβ0,
			c*cβ0-b*sβ0
	}
	D := math.Atan2(A[0], C[0])
	sD, cD := math.Sincos(D)
	for j := 1; j <= 8; j++ {
		X[j] = A[j]*cD - C[j]*sD
		Y[j] = A[j]*sD + C[j]*cD
		Z[j] = B[j]
		d := X[j] / s4[j].r
		X[j] += math.Abs(Z[j]) / k[j] * math.Sqrt(1-d*d)
		W := Δ / (Δ + Z[j]/2475)
		pos[j-1].X = X[j] * W
		pos[j-1].Y = Y[j] * W
	}
	return
}

// Physical computes quantities for physical observations of Jupiter.
//
// Results:
//	DS  Planetocentric declination of the Sun.
//	DE  Planetocentric declination of the Earth.
//	ω1  Longitude of the System I central meridian of the illuminated disk,
//	    as seen from Earth.
//	ω2  Longitude of the System II central meridian of the illuminated disk,
//	    as seen from Earth.
//	P   Geocentric position angle of Jupiter's northern rotation pole.
func Physical(jde float64, earth, jupiter *pp.V87Planet) (DS, DE, ω1, ω2, P unit.Angle) {
	// Step 1.
	d := jde - 2433282.5
	T1 := d / base.JulianCentury
	const p = math.Pi / 180
	α0 := 268*p + .1061*p*T1
	δ0 := 64.5*p - .0164*p*T1
	// Step 2.
	W1 := 17.71*p + 877.90003539*p*d
	W2 := 16.838*p + 870.27003539*p*d
	// Step 3.
	l0, b0, R := earth.Position(jde)
	l0, b0 = pp.ToFK5(l0, b0, jde)
	// Steps 4-7.
	sl0, cl0 := l0.Sincos()
	sb0 := b0.Sin()
	Δ := 4. // surely better than 0.
	var l, b unit.Angle
	var r, x, y, z float64
	f := func() {
		τ := base.LightTime(Δ)
		l, b, r = jupiter.Position(jde - τ)
		l, b = pp.ToFK5(l, b, jde)
		sb, cb := b.Sincos()
		sl, cl := l.Sincos()
		// (42.2) p. 289
		x = r*cb*cl - R*cl0
		y = r*cb*sl - R*sl0
		z = r*sb - R*sb0
		// (42.3) p. 289
		Δ = math.Sqrt(x*x + y*y + z*z)
	}
	f()
	f()
	// Step 8.
	ε0 := nutation.MeanObliquity(jde)
	// Step 9.
	sε0, cε0 := ε0.Sincos()
	sl, cl := l.Sincos()
	sb, cb := b.Sincos()
	αs := math.Atan2(cε0*sl-sε0*sb/cb, cl)
	δs := math.Asin(cε0*sb + sε0*cb*sl)
	// Step 10.
	sδs, cδs := math.Sincos(δs)
	sδ0, cδ0 := math.Sincos(δ0)
	DS = unit.Angle(math.Asin(-sδ0*sδs - cδ0*cδs*math.Cos(α0-αs)))
	// Step 11.
	u := y*cε0 - z*sε0
	v := y*sε0 + z*cε0
	α := math.Atan2(u, x)
	δ := math.Atan(v / math.Hypot(x, u))
	sδ, cδ := math.Sincos(δ)
	sα0α, cα0α := math.Sincos(α0 - α)
	ζ := math.Atan2(sδ0*cδ*cα0α-sδ*cδ0, cδ*sα0α)
	// Step 12.
	DE = unit.Angle(math.Asin(-sδ0*sδ - cδ0*cδ*math.Cos(α0-α)))
	// Step 13.
	ω1 = unit.Angle(W1 - ζ - 5.07033*p*Δ)
	ω2 = unit.Angle(W2 - ζ - 5.02626*p*Δ)
	// Step 14.
	C := unit.Angle((2*r*Δ + R*R - r*r - Δ*Δ) / (4 * r * Δ))
	if (l - l0).Sin() < 0 {
		C = -C
	}
	ω1 = (ω1 + C).Mod1()
	ω2 = (ω2 + C).Mod1()
	// Step 15.
	Δψ, Δε := nutation.Nutation(jde)
	ε := ε0 + Δε
	// Step 16.
	sε, cε := ε.Sincos()
	sα, cα := math.Sincos(α)
	α += .005693 * p * (cα*cl0*cε + sα*sl0) / cδ
	δ += .005693 * p * (cl0*cε*(sε/cε*cδ-sα*sδ) + cα*sδ*sl0)
	// Step 17.
	tδ := sδ / cδ
	Δα := (cε+sε*sα*tδ)*Δψ.Rad() - cα*tδ*Δε.Rad()
	Δδ := sε*cα*Δψ.Rad() + sα*Δε.Rad()
	αʹ := α + Δα
	δʹ := δ + Δδ
	sα0, cα0 := math.Sincos(α0)
	tδ0 := sδ0 / cδ0
	Δα0 := (cε+sε*sα0*tδ0)*Δψ.Rad() - cα0*tδ0*Δε.Rad()
	Δδ0 := sε*cα0*Δψ.Rad() + sα0*Δε.Rad()
	α0ʹ := α0 + Δα0
	δ0ʹ := δ0 + Δδ0
	// Step 18.
	sδʹ, cδʹ := math.Sincos(δʹ)
	sδ0ʹ, cδ0ʹ := math.Sincos(δ0ʹ)
	sα0ʹαʹ, cα0ʹαʹ := math.Sincos(α0ʹ - αʹ)
//.........这里部分代码省略.........

// Physical computes quantities for physical observations of Mars.
//
// Results:
//	DE  planetocentric declination of the Earth.
//	DS  planetocentric declination of the Sun.
//	ω   Areographic longitude of the central meridian, as seen from Earth.
//	P   Geocentric position angle of Mars' northern rotation pole.
//	Q   Position angle of greatest defect of illumination.
//	d   Apparent diameter of Mars.
//	q   Greatest defect of illumination.
//	k   Illuminated fraction of the disk.
func Physical(jde float64, earth, mars *pp.V87Planet) (DE, DS, ω, P, Q, d, q unit.Angle, k float64) {
	// Step 1.
	T := base.J2000Century(jde)
	const p = math.Pi / 180
	// (42.1) p. 288
	λ0 := 352.9065*p + 1.1733*p*T
	β0 := 63.2818*p - .00394*p*T
	// Step 2.
	l0, b0, R := earth.Position(jde)
	l0, b0 = pp.ToFK5(l0, b0, jde)
	// Steps 3, 4.
	sl0, cl0 := l0.Sincos()
	sb0 := b0.Sin()
	Δ := .5 // surely better than 0.
	τ := base.LightTime(Δ)
	var l, b unit.Angle
	var r, x, y, z float64
	f := func() {
		l, b, r = mars.Position(jde - τ)
		l, b = pp.ToFK5(l, b, jde)
		sb, cb := b.Sincos()
		sl, cl := l.Sincos()
		// (42.2) p. 289
		x = r*cb*cl - R*cl0
		y = r*cb*sl - R*sl0
		z = r*sb - R*sb0
		// (42.3) p. 289
		Δ = math.Sqrt(x*x + y*y + z*z)
		τ = base.LightTime(Δ)
	}
	f()
	f()
	// Step 5.
	λ := math.Atan2(y, x)
	β := math.Atan(z / math.Hypot(x, y))
	// Step 6.
	sβ0, cβ0 := math.Sincos(β0)
	sβ, cβ := math.Sincos(β)
	DE = unit.Angle(math.Asin(-sβ0*sβ - cβ0*cβ*math.Cos(λ0-λ)))
	// Step 7.
	N := 49.5581*p + .7721*p*T
	lʹ := l.Rad() - .00697*p/r
	bʹ := b.Rad() - .000225*p*math.Cos(l.Rad()-N)/r
	// Step 8.
	sbʹ, cbʹ := math.Sincos(bʹ)
	DS = unit.Angle(math.Asin(-sβ0*sbʹ - cβ0*cbʹ*math.Cos(λ0-lʹ)))
	// Step 9.
	W := 11.504*p + 350.89200025*p*(jde-τ-2433282.5)
	// Step 10.
	ε0 := nutation.MeanObliquity(jde)
	sε0, cε0 := ε0.Sincos()
	α0, δ0 := coord.EclToEq(unit.Angle(λ0), unit.Angle(β0), sε0, cε0)
	// Step 11.
	u := y*cε0 - z*sε0
	v := y*sε0 + z*cε0
	α := math.Atan2(u, x)
	δ := math.Atan(v / math.Hypot(x, u))
	sδ, cδ := math.Sincos(δ)
	sδ0, cδ0 := δ0.Sincos()
	sα0α, cα0α := math.Sincos(α0.Rad() - α)
	ζ := math.Atan2(sδ0*cδ*cα0α-sδ*cδ0, cδ*sα0α)
	// Step 12.
	ω = unit.Angle(W - ζ).Mod1()
	// Step 13.
	Δψ, Δε := nutation.Nutation(jde)
	// Step 14.
	sl0λ, cl0λ := math.Sincos(l0.Rad() - λ)
	λ += .005693 * p * cl0λ / cβ
	β += .005693 * p * sl0λ * sβ
	// Step 15.
	λ0 += Δψ.Rad()
	λ += Δψ.Rad()
	ε := ε0 + Δε
	// Step 16.
	sε, cε := ε.Sincos()
	α0ʹ, δ0ʹ := coord.EclToEq(unit.Angle(λ0), unit.Angle(β0), sε, cε)
	αʹ, δʹ := coord.EclToEq(unit.Angle(λ), unit.Angle(β), sε, cε)
	// Step 17.
	sδ0ʹ, cδ0ʹ := δ0ʹ.Sincos()
	sδʹ, cδʹ := δʹ.Sincos()
	sα0ʹαʹ, cα0ʹαʹ := (α0ʹ - αʹ).Sincos()
	// (42.4) p. 290
	P = unit.Angle(math.Atan2(cδ0ʹ*sα0ʹαʹ, sδ0ʹ*cδʹ-cδ0ʹ*sδʹ*cα0ʹαʹ))
	if P < 0 {
		P += 2 * math.Pi
	}
	// Step 18.
	s := l0 + math.Pi
	ss, cs := s.Sincos()
//.........这里部分代码省略.........

// Positions computes positions of moons of Jupiter.
//
// High accuracy method based on theory "E5."  Results returned in
// argument pos, which must not be nil.  Returned coordinates in units
// of Jupiter radii.
func E5(jde float64, earth, jupiter *pp.V87Planet, pos *[4]XY) {
	// I'll interject that I don't trust the results of this function.
	// There is obviously a great chance of typographic errors.
	// My Y results for the test case of the example don't agree with
	// Meeus's well at all, but do agree with the results from the less
	// accurate method.  This would seem to indicate a typo in Meeus's
	// computer implementation.  On the other hand, while my X results
	// agree reasonably well with his, our X results for satellite III
	// don't agree well with the result from the less accurate method,
	// perhaps indicating a typo in the presented algorithm.

	// variables assigned in following block
	var λ0, β0, t float64
	Δ := 5.
	{
		s, β, R := solar.TrueVSOP87(earth, jde)
		ss, cs := math.Sincos(s.Rad())
		sβ := math.Sin(β.Rad())
		τ := base.LightTime(Δ)
		var x, y, z float64
		f := func() {
			l, b, r := jupiter.Position(jde - τ)
			sl, cl := math.Sincos(l.Rad())
			sb, cb := math.Sincos(b.Rad())
			x = r*cb*cl + R*cs
			y = r*cb*sl + R*ss
			z = r*sb + R*sβ
			Δ = math.Sqrt(x*x + y*y + z*z)
			τ = base.LightTime(Δ)
		}
		f()
		f()
		λ0 = math.Atan2(y, x)
		β0 = math.Atan(z / math.Hypot(x, y))
		t = jde - 2443000.5 - τ
	}
	const p = math.Pi / 180
	l1 := 106.07719*p + 203.48895579*p*t
	l2 := 175.73161*p + 101.374724735*p*t
	l3 := 120.55883*p + 50.317609207*p*t
	l4 := 84.44459*p + 21.571071177*p*t
	π1 := 97.0881*p + .16138586*p*t
	π2 := 154.8663*p + .04726307*p*t
	π3 := 188.184*p + .00712734*p*t
	π4 := 335.2868*p + .00184*p*t
	ω1 := 312.3346*p - .13279386*p*t
	ω2 := 100.4411*p - .03263064*p*t
	ω3 := 119.1942*p - .00717703*p*t
	ω4 := 322.6186*p - .00175934*p*t
	Γ := .33033*p*math.Sin(163.679*p+.0010512*p*t) +
		.03439*p*math.Sin(34.486*p-.0161731*p*t)
	Φλ := 199.6766*p + .1737919*p*t
	ψ := 316.5182*p - .00000208*p*t
	G := 30.23756*p + .0830925701*p*t + Γ
	Gʹ := 31.97853*p + .0334597339*p*t
	const Π = 13.469942 * p

	Σ1 := .47259*p*math.Sin(2*(l1-l2)) +
		-.03478*p*math.Sin(π3-π4) +
		.01081*p*math.Sin(l2-2*l3+π3) +
		.00738*p*math.Sin(Φλ) +
		.00713*p*math.Sin(l2-2*l3+π2) +
		-.00674*p*math.Sin(π1+π3-2*Π-2*G) +
		.00666*p*math.Sin(l2-2*l3+π4) +
		.00445*p*math.Sin(l1-π3) +
		-.00354*p*math.Sin(l1-l2) +
		-.00317*p*math.Sin(2*ψ-2*Π) +
		.00265*p*math.Sin(l1-π4) +
		-.00186*p*math.Sin(G) +
		.00162*p*math.Sin(π2-π3) +
		.00158*p*math.Sin(4*(l1-l2)) +
		-.00155*p*math.Sin(l1-l3) +
		-.00138*p*math.Sin(ψ+ω3-2*Π-2*G) +
		-.00115*p*math.Sin(2*(l1-2*l2+ω2)) +
		.00089*p*math.Sin(π2-π4) +
		.00085*p*math.Sin(l1+π3-2*Π-2*G) +
		.00083*p*math.Sin(ω2-ω3) +
		.00053*p*math.Sin(ψ-ω2)
	Σ2 := 1.06476*p*math.Sin(2*(l2-l3)) +
		.04256*p*math.Sin(l1-2*l2+π3) +
		.03581*p*math.Sin(l2-π3) +
		.02395*p*math.Sin(l1-2*l2+π4) +
		.01984*p*math.Sin(l2-π4) +
		-.01778*p*math.Sin(Φλ) +
		.01654*p*math.Sin(l2-π2) +
		.01334*p*math.Sin(l2-2*l3+π2) +
		.01294*p*math.Sin(π3-π4) +
		-.01142*p*math.Sin(l2-l3) +
		-.01057*p*math.Sin(G) +
		-.00775*p*math.Sin(2*(ψ-Π)) +
		.00524*p*math.Sin(2*(l1-l2)) +
		-.0046*p*math.Sin(l1-l3) +
		.00316*p*math.Sin(ψ-2*G+ω3-2*Π) +
		-.00203*p*math.Sin(π1+π3-2*Π-2*G) +
		.00146*p*math.Sin(ψ-ω3) +
//.........这里部分代码省略.........

// LongitudeJ2000 returns geometric longitude referenced to equinox J2000.
func LongitudeJ2000(e *pp.V87Planet, jde float64) (l float64) {
	l, _, _ = e.Position2000(jde)
	return base.PMod(l+math.Pi-.09033/3600*math.Pi/180, 2*math.Pi)
}

// LongitudeJ2000 returns geometric longitude referenced to equinox J2000.
func LongitudeJ2000(e *pp.V87Planet, jde float64) (l unit.Angle) {
	l, _, _ = e.Position2000(jde)
	return (l + math.Pi - unit.AngleFromSec(.09033)).Mod1()
}