// EqToEcl converts equatorial coordinates to ecliptic coordinates.
func (ecl *Ecliptic) EqToEcl(eq *Equatorial, ε *Obliquity) *Ecliptic {
	sα, cα := eq.RA.Sincos()
	sδ, cδ := eq.Dec.Sincos()
	ecl.Lon = unit.Angle(math.Atan2(sα*ε.C+(sδ/cδ)*ε.S, cα)) // (13.1) p. 93
	ecl.Lat = unit.Angle(math.Asin(sδ*ε.C - cδ*ε.S*sα))      // (13.2) p. 93
	return ecl
}

// 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) (α unit.RA, δ unit.Angle) {
	L0, B0, R0 := earth.Position(jde)
	L, B, R := p.Position(jde)
	sB0, cB0 := B0.Sincos()
	sL0, cL0 := L0.Sincos()
	sB, cB := B.Sincos()
	sL, cL := L.Sincos()
	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 = B.Sincos()
		sL, cL = L.Sincos()
		x = R*cB*cL - R0*cB0*cL0
		y = R*cB*sL - R0*cB0*sL0
		z = R*sB - R0*sB0
	}
	λ := unit.Angle(math.Atan2(y, x))                // (33.1) p. 223
	β := unit.Angle(math.Atan2(z, math.Hypot(x, y))) // (33.2) p. 223
	Δλ, Δβ := apparent.EclipticAberration(λ, β, jde)
	λ, β = pp.ToFK5(λ+Δλ, β+Δβ, jde)
	Δψ, Δε := nutation.Nutation(jde)
	λ += Δψ
	sε, cε := (nutation.MeanObliquity(jde) + Δε).Sincos()
	return coord.EclToEq(λ, β, sε, cε)
	// Meeus gives a formula for elongation but doesn't spell out how to
	// obtain term λ0 and doesn't give an example solution.
}

// AngleError returns both an angle as in the function Angle, and an error
// as in the function Error.
//
// The algorithm is by B. Pessens.
func AngleError(r1, d1, r2, d2, r3, d3 unit.Angle) (ψ, ω unit.Angle) {
	sr1, cr1 := r1.Sincos()
	sd1, cd1 := d1.Sincos()
	sr2, cr2 := r2.Sincos()
	sd2, cd2 := d2.Sincos()
	sr3, cr3 := r3.Sincos()
	sd3, cd3 := d3.Sincos()
	a1 := cd1 * cr1
	a2 := cd2 * cr2
	a3 := cd3 * cr3
	b1 := cd1 * sr1
	b2 := cd2 * sr2
	b3 := cd3 * sr3
	c1 := sd1
	c2 := sd2
	c3 := sd3
	l1 := b1*c2 - b2*c1
	l2 := b2*c3 - b3*c2
	l3 := b1*c3 - b3*c1
	m1 := c1*a2 - c2*a1
	m2 := c2*a3 - c3*a2
	m3 := c1*a3 - c3*a1
	n1 := a1*b2 - a2*b1
	n2 := a2*b3 - a3*b2
	n3 := a1*b3 - a3*b1
	ψ = unit.Angle(math.Acos((l1*l2 + m1*m2 + n1*n2) /
		(math.Sqrt(l1*l1+m1*m1+n1*n1) * math.Sqrt(l2*l2+m2*m2+n2*n2))))
	ω = unit.Angle(math.Asin((a2*l3 + b2*m3 + c2*n3) /
		(math.Sqrt(a2*a2+b2*b2+c2*c2) * math.Sqrt(l3*l3+m3*m3+n3*n3))))
	return
}

// NewEclipticPrecessor constructs an EclipticPrecessor object and initializes
// it to precess coordinates from epochFrom to epochTo.
func NewEclipticPrecessor(epochFrom, epochTo float64) *EclipticPrecessor {
	// (21.5) p. 136
	ηCoeff := ηt
	πCoeff := πt
	pCoeff := pt
	if epochFrom != 2000 {
		T := (epochFrom - 2000) * .01
		ηCoeff = []float64{
			base.Horner(T, ηT...),
			-0.03302*s + 0.000598*s*T,
			0.000060 * s}
		πCoeff = []float64{
			base.Horner(T, πT...),
			-869.8089*s - 0.50491*s*T,
			0.03536 * s}
		pCoeff = []float64{
			base.Horner(T, pT...),
			1.11113*s - 0.000042*s*T,
			-0.000006 * s}
	}
	t := (epochTo - epochFrom) * .01
	p := &EclipticPrecessor{
		π: unit.Angle(base.Horner(t, πCoeff...)),
		p: unit.Angle(base.Horner(t, pCoeff...) * t),
	}
	η := unit.Angle(base.Horner(t, ηCoeff...) * t)
	p.sη, p.cη = η.Sincos()
	return p
}

// EqToHz computes Horizontal coordinates from equatorial coordinates.
//
//	α: right ascension coordinate to transform
//	δ: declination coordinate to transform
//	φ: latitude of observer on Earth
//	ψ: longitude of observer on Earth
//	st: sidereal time at Greenwich at time of observation.
//
// Sidereal time must be consistent with the equatorial coordinates.
// If coordinates are apparent, sidereal time must be apparent as well.
//
// Results:
//
//	A: azimuth of observed point, measured westward from the South.
//	h: elevation, or height of observed point above horizon.
func EqToHz(α unit.RA, δ, φ, ψ unit.Angle, st unit.Time) (A, h unit.Angle) {
	H := st.Rad() - ψ.Rad() - α.Rad()
	sH, cH := math.Sincos(H)
	sφ, cφ := φ.Sincos()
	sδ, cδ := ψ.Sincos()
	A = unit.Angle(math.Atan2(sH, cH*sφ-(sδ/cδ)*cφ)) // (13.5) p. 93
	h = unit.Angle(math.Asin(sφ*sδ + cφ*cδ*cH))      // (13.6) p. 93
	return
}

// Sep returns the angular separation between two celestial bodies.
//
// The algorithm is numerically naïve, and while patched up a bit for
// small separations, remains unstable for separations near π.
func Sep(r1, d1, r2, d2 unit.Angle) unit.Angle {
	sd1, cd1 := d1.Sincos()
	sd2, cd2 := d2.Sincos()
	cd := sd1*sd2 + cd1*cd2*(r1-r2).Cos() // (17.1) p. 109
	if cd < base.CosSmallAngle {
		return unit.Angle(math.Acos(cd))
	}
	// (17.2) p. 109
	return unit.Angle(math.Hypot((r2-r1).Rad()*cd1, (d2 - d1).Rad()))
}

func (m *moon) optical(λ, β unit.Angle) (lʹ, bʹ, A unit.Angle) {
	// (53.1) p. 372
	W := λ - m.Ω // (λ without nutation)
	sW, cW := W.Sincos()
	sβ, cβ := β.Sincos()
	A = unit.Angle(math.Atan2(sW*cβ*cI-sβ*sI, cW*cβ))
	lʹ = (A - m.F).Mod1()
	bʹ = unit.Angle(math.Asin(-sW*cβ*sI - sβ*cI))
	return
}

// Position computes apparent position angle and angular distance of
// components of a binary star.
//
//	e is eccentricity of the true orbit
//	a is angular apparent semimajor axis
//	i is inclination relative to the line of sight
//	Ω is position angle of the ascending node
//	ω is longitude of periastron
//	E is eccentric anomaly, computed for example with package kepler
//	   and the mean anomaly as returned by function M in this package.
//
// Return value θ is the apparent position angle, ρ is the angular distance.
func Position(e float64, a, i, Ω, ω, E unit.Angle) (θ, ρ unit.Angle) {
	r := a.Mul(1 - e*E.Cos())
	ν := unit.Angle(2 * math.Atan(math.Sqrt((1+e)/(1-e))*E.Div(2).Tan()))
	sνω, cνω := (ν + ω).Sincos()
	ci := i.Cos()
	num := sνω * ci
	θ = (unit.Angle(math.Atan2(num, cνω)) + Ω).Mod1()
	ρ = r.Mul(math.Sqrt(num*num + cνω*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
}

// EqToGal converts equatorial coordinates to galactic coordinates.
//
// Equatorial coordinates must be referred to the standard equinox of B1950.0.
// For conversion to B1950, see package precess and utility functions in
// package "unit".
func (g *Galactic) EqToGal(eq *Equatorial) *Galactic {
	sdα, cdα := (galacticNorth.RA - eq.RA).Sincos()
	sgδ, cgδ := galacticNorth.Dec.Sincos()
	sδ, cδ := eq.Dec.Sincos()
	// (13.7) p. 94
	x := unit.Angle(math.Atan2(sdα, cdα*sgδ-(sδ/cδ)*cgδ))
	g.Lon = (galacticLon0 + math.Pi - x).Mod1()
	// (13.8) p. 94
	g.Lat = unit.Angle(math.Asin(sδ*sgδ + cδ*cgδ*cdα))
	return g
}

// EclipticAtHorizon computes how the plane of the ecliptic intersects
// the horizon at a given local sidereal time as observed from a given
// geographic latitude.
//
//	ε is obliquity of the ecliptic.
//	φ is geographic latitude of observer.
//	θ is local sidereal time.
//
//	λ1 and λ2 are ecliptic longitudes where the ecliptic intersects the horizon.
//	I is the angle at which the ecliptic intersects the horizon.
func EclipticAtHorizon(ε, φ unit.Angle, θ unit.Time) (λ1, λ2, I unit.Angle) {
	sε, cε := ε.Sincos()
	sφ, cφ := φ.Sincos()
	sθ, cθ := θ.Angle().Sincos()
	// (14.2) p. 99
	λ := unit.Angle(math.Atan2(-cθ, sε*(sφ/cφ)+cε*sθ))
	if λ < 0 {
		λ += math.Pi
	}
	// (14.3) p. 99
	return λ, λ + math.Pi, unit.Angle(math.Acos(cε*sφ - sε*cφ*sθ))
}

// ReduceB1950ToJ2000 reduces orbital elements of a solar system body from
// equinox B1950 in the FK4 system to equinox J2000 in the FK5 system.
func ReduceB1950FK4ToJ2000FK5(eFrom, eTo *Elements) *Elements {
	W := _L + eFrom.Node
	si, ci := eFrom.Inc.Sincos()
	sJ, cJ := _J.Sincos()
	sW, cW := W.Sincos()
	eTo.Inc = unit.Angle(math.Acos(ci*cJ - si*sJ*cW))
	eTo.Node = (unit.Angle(math.Atan2(si*sW, ci*sJ+si*cJ*cW)) -
		_Lp).Mod1()
	eTo.Peri = (eFrom.Peri +
		unit.Angle(math.Atan2(sJ*sW, si*cJ+ci*sJ*cW))).Mod1()
	return eTo
}

// ReduceElements reduces orbital elements of a solar system body from one
// equinox to another.
//
// This function is described in chapter 24, but is located in this
// package so it can be a method of EclipticPrecessor.
func (p *EclipticPrecessor) ReduceElements(eFrom, eTo *elementequinox.Elements) *elementequinox.Elements {
	ψ := p.π + p.p
	si, ci := eFrom.Inc.Sincos()
	snp, cnp := (eFrom.Node - p.π).Sincos()
	// (24.1) p. 159
	eTo.Inc = unit.Angle(math.Acos(ci*p.cη + si*p.sη*cnp))
	// (24.2) p. 159
	eTo.Node = ψ +
		unit.Angle(math.Atan2(si*snp, p.cη*si*cnp-p.sη*ci))
	// (24.3) p. 160
	eTo.Peri = eFrom.Peri +
		unit.Angle(math.Atan2(-p.sη*snp, si*p.cη-ci*p.sη*cnp))
	return eTo
}

// Topocentric returns topocentric positions including parallax.
//
// Arguments α, δ are geocentric right ascension and declination in radians.
// Δ is distance to the observed object in AU.  ρsφʹ, ρcφʹ are parallax
// constants (see package globe.) L is geographic longitude of the observer,
// jde is time of observation.
//
// Results are observed topocentric ra and dec in radians.
func Topocentric(α unit.RA, δ unit.Angle, Δ, ρsφʹ, ρcφʹ float64, L unit.Angle, jde float64) (αʹ unit.RA, δʹ unit.Angle) {
	π := Horizontal(Δ)
	θ0 := sidereal.Apparent(jde)
	H := (θ0.Angle() - L - unit.Angle(α)).Mod1()
	sπ := π.Sin()
	sH, cH := H.Sincos()
	sδ, cδ := δ.Sincos()
	// (40.2) p. 279
	Δα := unit.HourAngle(math.Atan2(-ρcφʹ*sπ*sH, cδ-ρcφʹ*sπ*cH))
	αʹ = α.Add(Δα)
	// (40.3) p. 279
	δʹ = unit.Angle(math.Atan2((sδ-ρsφʹ*sπ)*Δα.Cos(), cδ-ρcφʹ*sπ*cH))
	return
}

// EclipticPrecess precesses coordinates eclFrom, leaving result in eclTo.
//
// The same struct may be used for eclFrom and eclTo.
// EclTo is returned for convenience.
func (p *EclipticPrecessor) Precess(eclFrom, eclTo *coord.Ecliptic) *coord.Ecliptic {
	// (21.7) p. 137
	sβ, cβ := eclFrom.Lat.Sincos()
	sd, cd := (p.π - eclFrom.Lon).Sincos()
	A := p.cη*cβ*sd - p.sη*sβ
	B := cβ * cd
	C := p.cη*sβ + p.sη*cβ*sd
	eclTo.Lon = p.p + p.π - unit.Angle(math.Atan2(A, B))
	if C < base.CosSmallAngle {
		eclTo.Lat = unit.Angle(math.Asin(C))
	} else {
		eclTo.Lat = unit.Angle(math.Acos(math.Hypot(A, B))) // near pole
	}
	return eclTo
}

// Topocentric3 returns topocentric hour angle and declination including parallax.
//
// This function implements the "alternative" method described in the text.
// The method should be similarly rigorous to that of Topocentric() and results
// should be virtually consistent.
func Topocentric3(α unit.RA, δ unit.Angle, Δ, ρsφʹ, ρcφʹ float64, L unit.Angle, jde float64) (Hʹ unit.HourAngle, δʹ unit.Angle) {
	π := Horizontal(Δ)
	θ0 := sidereal.Apparent(jde)
	H := (θ0.Angle() - L - unit.Angle(α)).Mod1()
	sπ := π.Sin()
	sH, cH := H.Sincos()
	sδ, cδ := δ.Sincos()
	A := cδ * sH
	B := cδ*cH - ρcφʹ*sπ
	C := sδ - ρsφʹ*sπ
	q := math.Sqrt(A*A + B*B + C*C)
	Hʹ = unit.HourAngle(math.Atan2(A, B))
	δʹ = unit.Angle(math.Asin(C / q))
	return
}

// EclToEq converts ecliptic coordinates to equatorial coordinates.
//
//	λ: ecliptic longitude coordinate to transform
//	β: ecliptic latitude coordinate to transform
//	sε: sine of obliquity of the ecliptic
//	cε: cosine of obliquity of the ecliptic
//
// Results:
//	α: right ascension
//	δ: declination
func EclToEq(λ, β unit.Angle, sε, cε float64) (α unit.RA, δ unit.Angle) {
	sλ, cλ := λ.Sincos()
	sβ, cβ := β.Sincos()
	α = unit.RAFromRad(math.Atan2(sλ*cε-(sβ/cβ)*sε, cλ)) // (13.3) p. 93
	δ = unit.Angle(math.Asin(sβ*cε + cβ*sε*sλ))          // (13.4) p. 93
	return
}

// AnomalyDistance returns true anomaly and distance for near-parabolic orbits.
//
// Distance r returned in AU.
// An error is returned if the algorithm fails to converge.
func (e *Elements) AnomalyDistance(jde float64) (ν unit.Angle, r float64, err error) {
	// fairly literal translation of code on p. 246
	q1 := base.K * math.Sqrt((1+e.Ecc)/e.PDis) / (2 * e.PDis) // line 20
	g := (1 - e.Ecc) / (1 + e.Ecc)                            // line 20

	t := jde - e.TimeP // line 22
	if t == 0 {        // line 24
		return 0, e.PDis, nil
	}
	d1, d := 10000., 1e-9        // line 14
	q2 := q1 * t                 // line 28
	s := 2. / (3 * math.Abs(q2)) // line 30
	s = 2 / math.Tan(2*math.Atan(math.Cbrt(math.Tan(math.Atan(s)/2))))
	if t < 0 { // line 34
		s = -s
	}
	if e.Ecc != 1 { // line 36
		l := 0 // line 38
		for {
			s0 := s // line 40
			z := 1.
			y := s * s
			g1 := -y * s
			q3 := q2 + 2*g*s*y/3 // line 42
			for {
				z += 1                          // line 44
				g1 = -g1 * g * y                // line 46
				z1 := (z - (z+1)*g) / (2*z + 1) // line 48
				f := z1 * g1                    // line 50
				q3 += f                         // line 52
				if z > 50 || math.Abs(f) > d1 { // line 54
					return 0, 0, errors.New("No convergence")
				}
				if math.Abs(f) <= d { // line 56
					break
				}
			}
			l++ // line 58
			if l > 50 {
				return 0, 0, errors.New("No convergence")
			}
			for {
				s1 := s // line 60
				s = (2*s*s*s/3 + q3) / (s*s + 1)
				if math.Abs(s-s1) <= d { // line 62
					break
				}
			}
			if math.Abs(s-s0) <= d { // line 64
				break
			}
		}
	}
	ν = unit.Angle(2 * math.Atan(s))               // line 66
	r = e.PDis * (1 + e.Ecc) / (1 + e.Ecc*ν.Cos()) // line 68
	if ν < 0 {                                     // line 70
		ν += 2 * math.Pi
	}
	return
}

// ReduceB1950ToJ2000 reduces orbital elements of a solar system body from
// equinox B1950 to J2000.
func ReduceB1950ToJ2000(eFrom, eTo *Elements) *Elements {
	// (24.4) p. 161
	const S = .0001139788
	const C = .9999999935
	W := eFrom.Node - unit.AngleFromDeg(174.298782)
	si, ci := eFrom.Inc.Sincos()
	sW, cW := W.Sincos()
	A := si * sW
	B := C*si*cW - S*ci
	eTo.Inc = unit.Angle(math.Asin(math.Hypot(A, B)))
	eTo.Node = (unit.AngleFromDeg(174.997194) +
		unit.Angle(math.Atan2(A, B))).Mod1()
	eTo.Peri = (eFrom.Peri +
		unit.Angle(math.Atan2(-S*sW, C*si-S*ci*cW))).Mod1()
	return eTo
}

func (m *moon) sun(λ, β unit.Angle, Δ float64, earth *pp.V87Planet) (l0, b0 unit.Angle) {
	λ0, _, R := solar.ApparentVSOP87(earth, m.jde)
	ΔR := unit.Angle(Δ / (R * base.AU))
	λH := λ0 + math.Pi + ΔR.Mul(β.Cos()*(λ0-λ).Sin())
	βH := ΔR * β
	return m.lib(λH, βH)
}