// Section "Trigonometric functions of large angles":
//
// Meeus makes his point, but an example with integer values is a bit unfair
// when trigonometric functions inherently work on floating point numbers.
func ExamplePMod_integer() {
	const large = 36000030

	// The direct function call loses precition as expected.
	fmt.Println("Direct:    ", math.Sin(large*math.Pi/180))

	// When the value is manually reduced to the integer 30, the Go constant
	// evaluaton does a good job of delivering a radian value to math.Sin that
	// evaluates to .5 exactly.
	fmt.Println("Integer 30:", math.Sin(30*math.Pi/180))

	// Math.Mod takes float64s and returns float64s.  The integer constants
	// here however can be represented exactly as float64s, and the returned
	// result is exact as well.
	fmt.Println("Math.Mod:  ", math.Mod(large, 360))

	// But when math.Mod is substituted into the Sin function, float64s
	// are multiplied instead of the high precision constants, and the result
	// comes back slightly inexact.
	fmt.Println("Sin Mod:   ", math.Sin(math.Mod(large, 360)*math.Pi/180))

	// Use of PMod on integer constants produces results identical to above.
	fmt.Println("PMod int:  ", math.Sin(base.PMod(large, 360)*math.Pi/180))

	// As soon as the large integer is scaled to a non-integer value though,
	// precision is lost and PMod is of no help recovering at this point.
	fmt.Println("PMod float:", math.Sin(base.PMod(large*math.Pi/180, 2*math.Pi)))
	// Output:
	// Direct:     0.49999999995724154
	// Integer 30: 0.5
	// Math.Mod:   30
	// Sin Mod:    0.49999999999999994
	// PMod int:   0.49999999999999994
	// PMod float: 0.49999999997845307
}

// 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, ☉, in radians
//	ν = true anomaly in radians
func True(T float64) (s, ν float64) {
	// (25.2) p. 163
	L0 := base.Horner(T, 280.46646, 36000.76983, 0.0003032) *
		math.Pi / 180
	M := MeanAnomaly(T)
	C := (base.Horner(T, 1.914602, -0.004817, -.000014)*
		math.Sin(M) +
		(0.019993-.000101*T)*math.Sin(2*M) +
		0.000289*math.Sin(3*M)) * math.Pi / 180
	return base.PMod(L0+C, 2*math.Pi), base.PMod(M+C, 2*math.Pi)
}

// 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 - 174.298782*math.Pi/180
	si, ci := math.Sincos(eFrom.Inc)
	sW, cW := math.Sincos(W)
	A := si * sW
	B := C*si*cW - S*ci
	eTo.Inc = math.Asin(math.Hypot(A, B))
	eTo.Node = base.PMod(174.997194*math.Pi/180+math.Atan2(A, B),
		2*math.Pi)
	eTo.Peri = base.PMod(eFrom.Peri+math.Atan2(-S*sW, C*si-S*ci*cW),
		2*math.Pi)
	return eTo
}

func TestNode0(t *testing.T) {
	for i, j := range n0 {
		if e := math.Abs(base.PMod(moon.Node(j)+1, 2*math.Pi) - 1); e > 1e-3 {
			t.Error(i, e)
		}
	}
}

// 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 {
	const (
		Lp = 4.50001688 * math.Pi / 180
		L  = 5.19856209 * math.Pi / 180
		J  = .00651966 * math.Pi / 180
	)
	W := L + eFrom.Node
	si, ci := math.Sincos(eFrom.Inc)
	sJ, cJ := math.Sincos(J)
	sW, cW := math.Sincos(W)
	eTo.Inc = math.Acos(ci*cJ - si*sJ*cW)
	eTo.Node = base.PMod(math.Atan2(si*sW, ci*sJ+si*cJ*cW)-Lp,
		2*math.Pi)
	eTo.Peri = base.PMod(eFrom.Peri+math.Atan2(sJ*sW, si*cJ+ci*sJ*cW),
		2*math.Pi)
	return eTo
}

// 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 = base.PMod(a.M0+k*a.M1, 360) * math.Pi / 180
	T = base.J2000Century(J)
	return
}

// Apparent0UT returns apparent sidereal time at Greenwich at 0h UT
// on the given JD.
//
// The result is in seconds of time and is in the range [0,86400).
func Apparent0UT(jd float64) float64 {
	j0, f := math.Modf(jd + .5)
	cen := (j0 - .5 - base.J2000) / 36525
	s := base.Horner(cen, iau82...) + f*1.00273790935*86400
	n := nutation.NutationInRA(j0)      // angle (radians) of RA
	ns := n * 3600 * 180 / math.Pi / 15 // convert RA to time in seconds
	return base.PMod(s+ns, 86400)
}

// GalToEq converts galactic coordinates to equatorial coordinates.
//
// Resulting equatorial coordinates will be referred to the standard equinox of
// B1950.0.  For subsequent conversion to other epochs, see package precess and
// utility functions in package meeus.
func GalToEq(l, b float64) (α, δ float64) {
	sdLon, cdLon := math.Sincos(l - galacticLon0)
	sgδ, cgδ := math.Sincos(galacticNorth.Dec)
	sb, cb := math.Sincos(b)
	y := math.Atan2(sdLon, cdLon*sgδ-(sb/cb)*cgδ)
	α = base.PMod(y+galacticNorth.RA, 2*math.Pi)
	δ = math.Asin(sb*sgδ + cb*cgδ*cdLon)
	return
}

// HzToEq transforms horizontal coordinates to equatorial coordinates.
//
//	A: azimuth
//	h: elevation
//	φ: 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:
//
//	α: right ascension
//	δ: declination
func HzToEq(A, h, φ, ψ, st float64) (α, δ float64) {
	sA, cA := math.Sincos(A)
	sh, ch := math.Sincos(h)
	sφ, cφ := math.Sincos(φ)
	H := math.Atan2(sA, cA*sφ+sh/ch*cφ)
	α = base.PMod(sexa.Time(st).Rad()-ψ-H, 2*math.Pi)
	δ = math.Asin(sφ*sh - cφ*ch*cA)
	return
}

// Topocentric2 returns topocentric corrections including parallax.
//
// This function implements the "non-rigorous" method descripted in the text.
//
// Note that results are corrections, not corrected coordinates.
func Topocentric2(α, δ, Δ, ρsφʹ, ρcφʹ, L, jde float64) (Δα, Δδ float64) {
	π := Horizontal(Δ)
	θ0 := base.Time(sidereal.Apparent(jde)).Rad()
	H := base.PMod(θ0-L-α, 2*math.Pi)
	sH, cH := math.Sincos(H)
	sδ, cδ := math.Sincos(δ)
	Δα = -π * ρcφʹ * sH / cδ         // (40.4) p. 280
	Δδ = -π * (ρsφʹ*cδ - ρcφʹ*cH*sδ) // (40.5) p. 280
	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 = 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
}

func (m *moon) optical(λ, β float64) (lʹ, bʹ, A float64) {
	// (53.1) p. 372
	W := λ - m.Ω // (λ without nutation)
	sW, cW := math.Sincos(W)
	sβ, cβ := math.Sincos(β)
	A = math.Atan2(sW*cβ*cI-sβ*sI, cW*cβ)
	lʹ = base.PMod(A-m.F, 2*math.Pi)
	bʹ = math.Asin(-sW*cβ*sI - sβ*cI)
	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
}

// 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(α, δ, Δ, ρsφʹ, ρcφʹ, L, jde float64) (αʹ, δʹ float64) {
	π := Horizontal(Δ)
	θ0 := base.Time(sidereal.Apparent(jde)).Rad()
	H := base.PMod(θ0-L-α, 2*math.Pi)
	sπ := math.Sin(π)
	sH, cH := math.Sincos(H)
	sδ, cδ := math.Sincos(δ)
	Δα := math.Atan2(-ρcφʹ*sπ*sH, cδ-ρcφʹ*sπ*cH) // (40.2) p. 279
	αʹ = α + Δα
	δʹ = math.Atan2((sδ-ρsφʹ*sπ)*math.Cos(Δα), cδ-ρcφʹ*sπ*cH) // (40.3) p. 279
	return
}

// Times computes UT rise, transit and set times for a celestial object on
// a day of interest.
//
// The function argurments do not actually include the day, but do include
// a number of values computed from the day.
//
//	p is geographic coordinates of observer.
//	ΔT is delta T.
//	h0 is "standard altitude" of the body.
//	Th0 is apparent sidereal time at 0h UT at Greenwich.
//	α3, δ3 are slices of three right ascensions and declinations.
//
// ΔT unit is seconds.  See package deltat.
//
// h0 unit is radians.
//
// Th0 must be the time on the day of interest, in seconds.
// See sidereal.Apparent0UT.
//
// α3, δ3 must be values at 0h dynamical time for the day before, the day of,
// and the day after the day of interest.  Units are radians.
//
// Result units are seconds and are in the range [0,86400)
func Times(p globe.Coord, ΔT, h0, Th0 float64, α3, δ3 []float64) (mRise, mTransit, mSet float64, err error) {
	mRise, mTransit, mSet, err = ApproxTimes(p, h0, Th0, α3[1], δ3[1])
	if err != nil {
		return
	}
	var d3α, d3δ *interp.Len3
	d3α, err = interp.NewLen3(-86400, 86400, α3)
	if err != nil {
		return
	}
	d3δ, err = interp.NewLen3(-86400, 86400, δ3)
	if err != nil {
		return
	}
	// adjust mTransit
	{
		th0 := base.PMod(Th0+mTransit*360.985647/360, 86400)
		α := d3α.InterpolateX(mTransit + ΔT)
		H := th0 - (p.Lon+α)*43200/math.Pi
		mTransit -= H
	}
	// adjust mRise, mSet
	sLat, cLat := math.Sincos(p.Lat)
	adjustRS := func(m float64) (float64, error) {
		th0 := base.PMod(Th0+m*360.985647/360, 86400)
		ut := m + ΔT
		α := d3α.InterpolateX(ut)
		δ := d3δ.InterpolateX(ut)
		H := th0 - (p.Lon+α)*43200/math.Pi
		sδ, cδ := math.Sincos(δ)
		h := sLat*sδ + cLat*cδ*math.Cos(H)
		return m + (h-h0)/cδ*cLat*math.Sin(H), nil
	}
	mRise, err = adjustRS(mRise)
	if err != nil {
		return
	}
	mSet, err = adjustRS(mSet)
	return
}

// ApproxTimes computes approximate UT rise, transit and set times for
// a celestial object on a day of interest.
//
// The function argurments do not actually include the day, but do include
// values computed from the day.
//
//	p is geographic coordinates of observer.
//	h0 is "standard altitude" of the body.
//	Th0 is apparent sidereal time at 0h UT at Greenwich.
//	α, δ are right ascension and declination of the body.
//
// h0 unit is radians.
//
// Th0 must be the time on the day of interest, in seconds.
// See sidereal.Apparent0UT.
//
// α, δ must be values at 0h dynamical time for the day of interest.
// Units are radians.
//
// Result units are seconds and are in the range [0,86400)
func ApproxTimes(p globe.Coord, h0, Th0 float64, α, δ float64) (mRise, mTransit, mSet float64, err error) {
	// Meeus works in a crazy mix of units.
	// This function and Times work with seconds of time as much as possible.

	// approximate local hour angle
	sLat, cLat := math.Sincos(p.Lat)
	sδ1, cδ1 := math.Sincos(δ)
	cH0 := (math.Sin(h0) - sLat*sδ1) / (cLat * cδ1) // (15.1) p. 102
	if cH0 < -1 || cH0 > 1 {
		err = ErrorCircumpolar
		return
	}
	H0 := math.Acos(cH0) * 43200 / math.Pi

	// approximate transit, rise, set times.
	// (15.2) p. 102.
	mt := (α+p.Lon)*43200/math.Pi - Th0
	mTransit = base.PMod(mt, 86400)
	mRise = base.PMod(mt-H0, 86400)
	mSet = base.PMod(mt+H0, 86400)
	return
}

// Section "Trigonometric functions of large angles":
//
// A non integer example better illustrates that reduction to a range
// does not rescue precision.
func ExamplePMod_mars() {
	// Angle W from step 9 of example 42.a, as suggested.
	const W = 5492522.4593
	const reduced = 2.4593

	// Direct function call.  It's a number.  How correct is it?
	fmt.Println("Direct:  ", math.Sin(W*math.Pi/180))

	// Manually reduced to range 0-360.  This is presumably the "correct"
	// answer, but note that the reduced number has a reduced number of
	// significat digits.  The answer cannot have any more significant digits.
	fmt.Println("Reduced: ", math.Sin(reduced*math.Pi/180))

	// Accordingly, PMod cannot rescue any precision, whether done on degrees
	// or radians.
	fmt.Println("PMod deg:", math.Sin(base.PMod(W, 360)*math.Pi/180))
	fmt.Println("PMod rad:", math.Sin(base.PMod(W*math.Pi/180, 2*math.Pi)))
	// Output:
	// Direct:   0.04290970350270464
	// Reduced:  0.04290970350923273
	// PMod deg: 0.04290970351307828
	// PMod rad: 0.04290970350643808
}

// Physical2 computes quantities for physical observations of Jupiter.
//
// Results are less accurate than with Physical().
//
// 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.
//
// All angular results in radians.
func Physical2(jde float64) (DS, DE, ω1, ω2 float64) {
	d := jde - base.J2000
	const p = math.Pi / 180
	V := 172.74*p + .00111588*p*d
	M := 357.529*p + .9856003*p*d
	sV := math.Sin(V)
	N := 20.02*p + .0830853*p*d + .329*p*sV
	J := 66.115*p + .9025179*p*d - .329*p*sV
	sM, cM := math.Sincos(M)
	sN, cN := math.Sincos(N)
	s2M, c2M := math.Sincos(2 * M)
	s2N, c2N := math.Sincos(2 * N)
	A := 1.915*p*sM + .02*p*s2M
	B := 5.555*p*sN + .168*p*s2N
	K := J + A - B
	R := 1.00014 - .01671*cM - .00014*c2M
	r := 5.20872 - .25208*cN - .00611*c2N
	sK, cK := math.Sincos(K)
	Δ := math.Sqrt(r*r + R*R - 2*r*R*cK)
	ψ := math.Asin(R / Δ * sK)
	dd := d - Δ/173
	ω1 = 210.98*p + 877.8169088*p*dd + ψ - B
	ω2 = 187.23*p + 870.1869088*p*dd + ψ - B
	C := math.Sin(ψ / 2)
	C *= C
	if sK > 0 {
		C = -C
	}
	ω1 = base.PMod(ω1+C, 2*math.Pi)
	ω2 = base.PMod(ω2+C, 2*math.Pi)
	λ := 34.35*p + .083091*p*d + .329*p*sV + B
	DS = 3.12 * p * math.Sin(λ+42.8*p)
	DE = DS - 2.22*p*math.Sin(ψ)*math.Cos(λ+22*p) -
		1.3*p*(r-Δ)/Δ*math.Sin(λ-100.5*p)
	return
}

// 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(α, δ, Δ, ρsφʹ, ρcφʹ, L, jde float64) (Hʹ, δʹ float64) {
	π := Horizontal(Δ)
	θ0 := base.Time(sidereal.Apparent(jde)).Rad()
	H := base.PMod(θ0-L-α, 2*math.Pi)
	sπ := math.Sin(π)
	sH, cH := math.Sincos(H)
	sδ, cδ := math.Sincos(δ)
	A := cδ * sH
	B := cδ*cH - ρcφʹ*sπ
	C := sδ - ρsφʹ*sπ
	q := math.Sqrt(A*A + B*B + C*C)
	Hʹ = math.Atan2(A, B)
	δʹ = math.Asin(C / q)
	return
}

// E computes the "equation of time" for the given JDE.
//
// Parameter e must be a planetposition.V87Planet object for Earth obtained
// with planetposition.LoadPlanet.
//
// Result is equation of time as an hour angle in radians.
func E(jde float64, e *pp.V87Planet) float64 {
	τ := base.J2000Century(jde) * .1
	L0 := l0(τ)
	// code duplicated from solar.ApparentEquatorialVSOP87 so that
	// we can keep Δψ and cε
	s, β, R := solar.TrueVSOP87(e, jde)
	Δψ, Δε := nutation.Nutation(jde)
	a := -20.4898 / 3600 * math.Pi / 180 / R
	λ := s + Δψ + a
	ε := nutation.MeanObliquity(jde) + Δε
	sε, cε := math.Sincos(ε)
	α, _ := coord.EclToEq(λ, β, sε, cε)
	// (28.1) p. 183
	E := L0 - .0057183*math.Pi/180 - α + Δψ*cε
	return base.PMod(E+math.Pi, 2*math.Pi) - math.Pi
}