Golang base.PMod函数代码示例

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