Golang Angle.Rad方法代码示例

// Kepler1 solves Kepler's equation by iteration.
//
// The iterated formula is
//
//	E1 = M + e * sin(E0)
//
// Argument e is eccentricity, M is mean anomaly,
// places is the desired number of decimal places in the result.
//
// Result E is eccentric anomaly.
//
// For some vaues of e and M it will fail to converge and the
// function will return an error.
func Kepler1(e float64, M unit.Angle, places int) (E unit.Angle, err error) {
	f := func(E0 float64) float64 {
		return M.Rad() + e*math.Sin(E0) // (30.5) p. 195
	}
	ea, err := iterate.DecimalPlaces(f, M.Rad(), places, places*5)
	return unit.Angle(ea), err
}

// Kepler2 solves Kepler's equation by iteration.
//
// The iterated formula is
//
//	E1 = E0 + (M + e * sin(E0) - E0) / (1 - e * cos(E0))
//
// Argument e is eccentricity, M is mean anomaly,
// places is the desired number of decimal places in the result.
//
// Result E is eccentric anomaly.
//
// The function converges over a wider range of inputs than does Kepler1
// but it also fails to converge for some values of e and M.
func Kepler2(e float64, M unit.Angle, places int) (E unit.Angle, err error) {
	f := func(E0 float64) float64 {
		se, ce := math.Sincos(E0)
		return E0 + (M.Rad()+e*se-E0)/(1-e*ce) // (30.7) p. 199
	}
	ea, err := iterate.DecimalPlaces(f, M.Rad(), places, places)
	return unit.Angle(ea), err
}

// Kepler2a solves Kepler's equation by iteration.
//
// The iterated formula is the same as in Kepler2 but a limiting function
// avoids divergence.
//
// Argument e is eccentricity, M is mean anomaly,
// places is the desired number of decimal places in the result.
//
// Result E is eccentric anomaly.
func Kepler2a(e float64, M unit.Angle, places int) (E unit.Angle, err error) {
	f := func(E0 float64) float64 {
		se, ce := math.Sincos(E0)
		// method of Leingärtner, p. 205
		return E0 + math.Asin(math.Sin((M.Rad()+e*se-E0)/(1-e*ce)))
	}
	ea, err := iterate.DecimalPlaces(f, M.Rad(), places, places*5)
	return unit.Angle(ea), err
}

// Kepler2b solves Kepler's equation by iteration.
//
// The iterated formula is the same as in Kepler2 but a (different) limiting
// function avoids divergence.
//
// Argument e is eccentricity, M is mean anomaly,
// places is the desired number of decimal places in the result.
//
// Result E is eccentric anomaly.
func Kepler2b(e float64, M unit.Angle, places int) (E unit.Angle, err error) {
	f := func(E0 float64) float64 {
		se, ce := math.Sincos(E0)
		d := (M.Rad() + e*se - E0) / (1 - e*ce)
		// method of Steele, p. 205
		if d > .5 {
			d = .5
		} else if d < -.5 {
			d = -.5
		}
		return E0 + d
	}
	ea, err := iterate.DecimalPlaces(f, M.Rad(), places, places)
	return unit.Angle(ea), err
}

// 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()
//.........这里部分代码省略.........

// ApproxLinearDistance computes a distance across the surface of the Earth.
//
// Approximating the Earth as a sphere, the function takes a geocentric angular
// distance and returns the corresponding linear distance in Km.
func ApproxLinearDistance(d unit.Angle) float64 {
	return 6371 * d.Rad()
}