Golang base.Horner函数代码示例

func newMp(y, q float64) *mp {
	m := &mp{k: snap(y, q)}
	m.T = m.k * ck // (49.3) p. 350
	m.E = base.Horner(m.T, 1, -.002516, -.0000074)
	m.M = base.Horner(m.T, 2.5534*p, 29.1053567*p/ck,
		-.0000014*p, -.00000011*p)
	m.Mʹ = base.Horner(m.T, 201.5643*p, 385.81693528*p/ck,
		.0107582*p, .00001238*p, -.000000058*p)
	m.F = base.Horner(m.T, 160.7108*p, 390.67050284*p/ck,
		-.0016118*p, -.00000227*p, .000000011*p)
	m.Ω = base.Horner(m.T, 124.7746*p, -1.56375588*p/ck,
		.0020672*p, .00000215*p)
	m.A[0] = 299.7*p + .107408*p*m.k - .009173*m.T*m.T
	m.A[1] = 251.88*p + .016321*p*m.k
	m.A[2] = 251.83*p + 26.651886*p*m.k
	m.A[3] = 349.42*p + 36.412478*p*m.k
	m.A[4] = 84.66*p + 18.206239*p*m.k
	m.A[5] = 141.74*p + 53.303771*p*m.k
	m.A[6] = 207.17*p + 2.453732*p*m.k
	m.A[7] = 154.84*p + 7.30686*p*m.k
	m.A[8] = 34.52*p + 27.261239*p*m.k
	m.A[9] = 207.19*p + .121824*p*m.k
	m.A[10] = 291.34*p + 1.844379*p*m.k
	m.A[11] = 161.72*p + 24.198154*p*m.k
	m.A[12] = 239.56*p + 25.513099*p*m.k
	m.A[13] = 331.55*p + 3.592518*p*m.k
	return m
}

// Nutation returns nutation in longitude (Δψ) and nutation in obliquity (Δε)
// for a given JDE.
//
// JDE = UT + ΔT, see package deltat.
//
// Computation is by 1980 IAU theory, with terms < .0003″ neglected.
func Nutation(jde float64) (Δψ, Δε unit.Angle) {
	T := base.J2000Century(jde)
	D := base.Horner(T,
		297.85036, 445267.11148, -0.0019142, 1./189474) * math.Pi / 180
	M := base.Horner(T,
		357.52772, 35999.050340, -0.0001603, -1./300000) * math.Pi / 180
	N := base.Horner(T,
		134.96298, 477198.867398, 0.0086972, 1./5620) * math.Pi / 180
	F := base.Horner(T,
		93.27191, 483202.017538, -0.0036825, 1./327270) * math.Pi / 180
	Ω := base.Horner(T,
		125.04452, -1934.136261, 0.0020708, 1./450000) * math.Pi / 180
	// sum in reverse order to accumulate smaller terms first
	var Δψs, Δεs float64
	for i := len(table22A) - 1; i >= 0; i-- {
		row := table22A[i]
		arg := row.d*D + row.m*M + row.n*N + row.f*F + row.ω*Ω
		s, c := math.Sincos(arg)
		Δψs += s * (row.s0 + row.s1*T)
		Δεs += c * (row.c0 + row.c1*T)
	}
	Δψ = unit.AngleFromSec(Δψs * .0001)
	Δε = unit.AngleFromSec(Δεs * .0001)
	return
}

// Extremum returns the x and y values at the extremum.
//
// Results are restricted to the range of the table given to the constructor
// NewLen5.  (Meeus actually recommends restricting the range to one unit of
// the tabular interval, but that seems a little harsh.)
func (d *Len5) Extremum() (x, y float64, err error) {
	// (3.9) p. 29
	nCoeff := []float64{
		6*(d.b+d.c) - d.h - d.j,
		0,
		3 * (d.h + d.k),
		2 * d.k,
	}
	den := d.k - 12*d.f
	if den == 0 {
		return 0, 0, ErrorExtremumOutside
	}
	n0, ok := iterate(0, func(n0 float64) float64 {
		return base.Horner(n0, nCoeff...) / den
	})
	if !ok {
		return 0, 0, ErrorNoConverge
	}
	if n0 < -2 || n0 > 2 {
		return 0, 0, ErrorExtremumOutside
	}
	x = .5*d.xSum + .25*d.xDiff*n0
	y = base.Horner(n0, d.interpCoeff...)
	return x, y, nil
}

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

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

func newLa(y, h float64) *la {
	l := &la{k: snap(y, h)}
	l.T = l.k * ck // (50.3) p. 350
	l.D = base.Horner(l.T, 171.9179*p, 335.9106046*p/ck,
		-.0100383*p, -.00001156*p, .000000055*p)
	l.M = base.Horner(l.T, 347.3477*p, 27.1577721*p/ck,
		-.000813*p, -.000001*p)
	l.F = base.Horner(l.T, 316.6109*p, 364.5287911*p/ck,
		-.0125053*p, -.0000148*p)
	return l
}

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

func dmf(T float64) (D, M, Mʹ, F float64) {
	D = base.Horner(T, 297.8501921*p, 445267.1114034*p,
		-.0018819*p, p/545868, -p/113065000)
	M = base.Horner(T, 357.5291092*p, 35999.0502909*p,
		-.0001535*p, p/24490000)
	Mʹ = base.Horner(T, 134.9633964*p, 477198.8675055*p,
		.0087414*p, p/69699, -p/14712000)
	F = base.Horner(T, 93.272095*p, 483202.0175233*p,
		-.0036539*p, -p/3526000, p/863310000)
	return
}

// Sum computes a sum of periodic terms.
func sum(T, M float64, c [][]float64) float64 {
	j := base.Horner(T, c[0]...)
	mm := 0.
	for i := 1; i < len(c); i++ {
		mm += M
		smm, cmm := math.Sincos(mm)
		j += smm * base.Horner(T, c[i]...)
		i++
		j += cmm * base.Horner(T, c[i]...)
	}
	return j
}

// SumA computes the sum of periodic terms with "additional angles"
func sumA(T, M float64, c [][]float64, aa []caa) float64 {
	i := len(c) - 2*len(aa)
	j := sum(T, M, c[:i])
	for k := 0; k < len(aa); k++ {
		saa, caa := math.Sincos((aa[k].c + aa[k].f*T) * math.Pi / 180)
		j += saa * base.Horner(T, c[i]...)
		i++
		j += caa * base.Horner(T, c[i]...)
		i++
	}
	return j
}

func (q *qs) iapetus() (r r4) {
	L := 261.1582*d + 22.57697855*d*q.t4
	ϖʹ := 91.796*d + .562*d*q.t7
	ψ := 4.367*d - .195*d*q.t7
	θ := 146.819*d - 3.198*d*q.t7
	φ := 60.47*d + 1.521*d*q.t7
	Φ := 205.055*d - 2.091*d*q.t7
	eʹ := .028298 + .001156*q.t11
	ϖ0 := 352.91*d + 11.71*d*q.t11
	μ := 76.3852*d + 4.53795125*d*q.t10
	iʹ := base.Horner(q.t11, 18.4602*d, -.9518*d, -.072*d, .0054*d)
	Ωʹ := base.Horner(q.t11, 143.198*d, -3.919*d, .116*d, .008*d)
	l := μ - ϖ0
	g := ϖ0 - Ωʹ - ψ
	g1 := ϖ0 - Ωʹ - φ
	ls := q.W5 - ϖʹ
	gs := ϖʹ - θ
	lT := L - q.W4
	gT := q.W4 - Φ
	u1 := 2 * (l + g - ls - gs)
	u2 := l + g1 - lT - gT
	u3 := l + 2*(g-ls-gs)
	u4 := lT + gT - g1
	u5 := 2 * (ls + gs)
	sl, cl := math.Sincos(l)
	su1, cu1 := math.Sincos(u1)
	su2, cu2 := math.Sincos(u2)
	su3, cu3 := math.Sincos(u3)
	su4, cu4 := math.Sincos(u4)
	slu2, clu2 := math.Sincos(l + u2)
	sg1gT, cg1gT := math.Sincos(g1 - gT)
	su52g, cu52g := math.Sincos(u5 - 2*g)
	su5ψ, cu5ψ := math.Sincos(u5 + ψ)
	su2φ, cu2φ := math.Sincos(u2 + φ)
	s5, c5 := math.Sincos(l + g1 + lT + gT + φ)
	a := 58.935028 + .004638*cu1 + .058222*cu2
	e := eʹ - .0014097*cg1gT + .0003733*cu52g +
		.000118*cu3 + .0002408*cl + .0002849*clu2 + .000619*cu4
	w := .08077*d*sg1gT + .02139*d*su52g - .00676*d*su3 +
		.0138*d*sl + .01632*d*slu2 + .03547*d*su4
	p := ϖ0 + w/eʹ
	λʹ := μ - .04299*d*su2 - .00789*d*su1 - .06312*d*math.Sin(ls) -
		.00295*d*math.Sin(2*ls) - .02231*d*math.Sin(u5) + .0065*d*su5ψ
	i := iʹ + .04204*d*cu5ψ + .00235*d*c5 + .0036*d*cu2φ
	wʹ := .04204*d*su5ψ + .00235*d*s5 + .00358*d*su2φ
	Ω := Ωʹ + wʹ/math.Sin(iʹ)
	return q.subr(λʹ, p, e, a, Ω, i)
}

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

// MeanObliquity returns mean obliquity (ε₀) following the IAU 1980
// polynomial.
//
// Accuracy is 1″ over the range 1000 to 3000 years and 10″ over the range
// 0 to 4000 years.
func MeanObliquity(jde float64) unit.Angle {
	// (22.2) p. 147
	return unit.AngleFromSec(base.Horner(base.J2000Century(jde),
		unit.FromSexaSec(' ', 23, 26, 21.448),
		-46.815,
		-0.00059,
		0.001813))
}

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

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