Golang Angle.Sincos方法代码示例

// Vertical computes data for a vertical sundial.
//
// Argument φ is geographic latitude at which the sundial will be located.
// D is gnomonic declination, the azimuth of the perpendicular to the plane
// of the sundial, measured from the southern meridian towards the west.
// Argument a is the length of a straight stylus perpendicular to the plane
// of the sundial.
//
// Results consist of a set of lines, a center point, and u, the length of a
// polar stylus.  They are in units of a, the stylus length.
func Vertical(φ, D unit.Angle, a float64) (lines []Line, center Point, u float64) {
	sφ, cφ := φ.Sincos()
	tφ := sφ / cφ
	sD, cD := D.Sincos()
	for i := 0; i < 24; i++ {
		l := Line{Hour: i}
		H := float64(i-12) * 15 * math.Pi / 180
		aH := math.Abs(H)
		sH, cH := math.Sincos(H)
		for _, d := range m {
			tδ := math.Tan(d * math.Pi / 180)
			H0 := math.Acos(-tφ * tδ)
			if aH > H0 {
				continue // sun below horizon
			}
			Q := sD*sH + sφ*cD*cH - cφ*cD*tδ
			if Q < 0 {
				continue // sun below plane of sundial
			}
			x := a * (cD*sH - sφ*sD*cH + cφ*sD*tδ) / Q
			y := -a * (cφ*cH + sφ*tδ) / Q
			l.Points = append(l.Points, Point{x, y})
		}
		if len(l.Points) > 0 {
			lines = append(lines, l)
		}
	}
	center.X = -a * sD / cD
	center.Y = a * tφ / cD
	u = a / math.Abs(cφ*cD)
	return
}

func (m *moon) physical(A, bʹ unit.Angle) (lʺ, bʺ unit.Angle) {
	// (53.2) p. 373
	sA, cA := A.Sincos()
	lʺ = -m.τ + (m.ρ.Mul(cA) + m.σ.Mul(sA)).Mul(bʹ.Tan())
	bʺ = m.σ.Mul(cA) - m.ρ.Mul(sA)
	return
}

// Angle returns the angle between great circles defined by three points.
//
// Coordinates may be right ascensions and declinations or longitudes and
// latitudes.  If r1, d1, r2, d2 defines one line and r2, d2, r3, d3 defines
// another, the result is the angle between the two lines.
//
// Algorithm by Meeus.
func Angle(r1, d1, r2, d2, r3, d3 unit.Angle) unit.Angle {
	sd2, cd2 := d2.Sincos()
	sr21, cr21 := (r2 - r1).Sincos()
	sr32, cr32 := (r3 - r2).Sincos()
	C1 := math.Atan2(sr21, cd2*d1.Tan()-sd2*cr21)
	C2 := math.Atan2(sr32, cd2*d3.Tan()-sd2*cr32)
	return unit.Angle(C1 + C2)
}

// SepPauwels returns the angular separation between two celestial bodies.
//
// The algorithm is a numerically stable form of that used in Sep.
func SepPauwels(r1, d1, r2, d2 unit.Angle) unit.Angle {
	sd1, cd1 := d1.Sincos()
	sd2, cd2 := d2.Sincos()
	cdr := (r2 - r1).Cos()
	x := cd1*sd2 - sd1*cd2*cdr
	y := cd2 * (r2 - r1).Sin()
	z := sd1*sd2 + cd1*cd2*cdr
	return unit.Angle(math.Atan2(math.Hypot(x, y), z))
}

// 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 unit.Angle) (α unit.RA, δ unit.Angle) {
	sdLon, cdLon := (l - galacticLon0).Sincos()
	sgδ, cgδ := galacticNorth.Dec.Sincos()
	sb, cb := b.Sincos()
	y := math.Atan2(sdLon, cdLon*sgδ-(sb/cb)*cgδ)
	α = unit.RAFromRad(y + galacticNorth.RA.Rad())
	δ = unit.Angle(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
// in the sense that tf coordinates are apparent, sidereal time must be
// apparent as well.
//
// Results:
//
//	α: right ascension
//	δ: declination
func HzToEq(A, h, φ, ψ unit.Angle, st unit.Time) (α unit.RA, δ unit.Angle) {
	sA, cA := A.Sincos()
	sh, ch := h.Sincos()
	sφ, cφ := φ.Sincos()
	H := math.Atan2(sA, cA*sφ+sh/ch*cφ)
	α = unit.RAFromRad(st.Rad() - ψ.Rad() - H)
	δ = unit.Angle(math.Asin(sφ*sh - cφ*ch*cA))
	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()))
}

// General computes data for the general case of a planar sundial.
//
// Argument φ is geographic latitude at which the sundial will be located.
// D is gnomonic declination, the azimuth of the perpendicular to the plane
// of the sundial, measured from the southern meridian towards the west.
// Argument a is the length of a straight stylus perpendicular to the plane
// of the sundial, z is zenithal distance of the direction defined by the
// stylus.  Units of stylus length a are arbitrary.
//
// Results consist of a set of lines, a center point, u, the length of a
// polar stylus, and ψ, the angle which the polar stylus makes with the plane
// of the sundial.  The center point, the points defining the hour lines, and
// u are in units of a, the stylus length.
func General(φ, D unit.Angle, a float64, z unit.Angle) (lines []Line, center Point, u float64, ψ unit.Angle) {
	sφ, cφ := φ.Sincos()
	tφ := sφ / cφ
	sD, cD := D.Sincos()
	sz, cz := z.Sincos()
	P := sφ*cz - cφ*sz*cD
	for i := 0; i < 24; i++ {
		l := Line{Hour: i}
		H := float64(i-12) * 15 * math.Pi / 180
		aH := math.Abs(H)
		sH, cH := math.Sincos(H)
		for _, d := range m {
			tδ := math.Tan(d * math.Pi / 180)
			H0 := math.Acos(-tφ * tδ)
			if aH > H0 {
				continue // sun below horizon
			}
			Q := sD*sz*sH + (cφ*cz+sφ*sz*cD)*cH + P*tδ
			if Q < 0 {
				continue // sun below plane of sundial
			}
			Nx := cD*sH - sD*(sφ*cH-cφ*tδ)
			Ny := cz*sD*sH - (cφ*sz-sφ*cz*cD)*cH - (sφ*sz+cφ*cz*cD)*tδ
			l.Points = append(l.Points, Point{a * Nx / Q, a * Ny / Q})
		}
		if len(l.Points) > 0 {
			lines = append(lines, l)
		}
	}
	center.X = a / P * cφ * sD
	center.Y = -a / P * (sφ*sz + cφ*cz*cD)
	aP := math.Abs(P)
	u = a / aP
	ψ = unit.Angle(math.Asin(aP))
	return
}

// Error returns an error angle of three nearly co-linear points.
//
// For the line defined by r1, d1, r2, d2, the result is the anglular distance
// between that line and r0, d0.
//
// Algorithm by Meeus.
func Error(r1, d1, r2, d2, r0, d0 unit.Angle) unit.Angle {
	sr1, cr1 := r1.Sincos()
	sd1, cd1 := d1.Sincos()
	sr2, cr2 := r2.Sincos()
	sd2, cd2 := d2.Sincos()
	X1 := cd1 * cr1
	X2 := cd2 * cr2
	Y1 := cd1 * sr1
	Y2 := cd2 * sr2
	Z1 := sd1
	Z2 := sd2
	A := Y1*Z2 - Z1*Y2
	B := Z1*X2 - X1*Z2
	C := X1*Y2 - Y1*X2
	m := r0.Tan()
	n := d0.Tan() / r0.Cos()
	return unit.Angle(math.Asin((A + B*m + C*n) /
		(math.Sqrt(A*A+B*B+C*C) * math.Sqrt(1+m*m+n*n))))
}

// Physical computes quantities for physical observations of Jupiter.
//
// 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.
//	P   Geocentric position angle of Jupiter's northern rotation pole.
func Physical(jde float64, earth, jupiter *pp.V87Planet) (DS, DE, ω1, ω2, P unit.Angle) {
	// Step 1.
	d := jde - 2433282.5
	T1 := d / base.JulianCentury
	const p = math.Pi / 180
	α0 := 268*p + .1061*p*T1
	δ0 := 64.5*p - .0164*p*T1
	// Step 2.
	W1 := 17.71*p + 877.90003539*p*d
	W2 := 16.838*p + 870.27003539*p*d
	// Step 3.
	l0, b0, R := earth.Position(jde)
	l0, b0 = pp.ToFK5(l0, b0, jde)
	// Steps 4-7.
	sl0, cl0 := l0.Sincos()
	sb0 := b0.Sin()
	Δ := 4. // surely better than 0.
	var l, b unit.Angle
	var r, x, y, z float64
	f := func() {
		τ := base.LightTime(Δ)
		l, b, r = jupiter.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)
	}
	f()
	f()
	// Step 8.
	ε0 := nutation.MeanObliquity(jde)
	// Step 9.
	sε0, cε0 := ε0.Sincos()
	sl, cl := l.Sincos()
	sb, cb := b.Sincos()
	αs := math.Atan2(cε0*sl-sε0*sb/cb, cl)
	δs := math.Asin(cε0*sb + sε0*cb*sl)
	// Step 10.
	sδs, cδs := math.Sincos(δs)
	sδ0, cδ0 := math.Sincos(δ0)
	DS = unit.Angle(math.Asin(-sδ0*sδs - cδ0*cδs*math.Cos(α0-αs)))
	// 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α := math.Sincos(α0 - α)
	ζ := math.Atan2(sδ0*cδ*cα0α-sδ*cδ0, cδ*sα0α)
	// Step 12.
	DE = unit.Angle(math.Asin(-sδ0*sδ - cδ0*cδ*math.Cos(α0-α)))
	// Step 13.
	ω1 = unit.Angle(W1 - ζ - 5.07033*p*Δ)
	ω2 = unit.Angle(W2 - ζ - 5.02626*p*Δ)
	// Step 14.
	C := unit.Angle((2*r*Δ + R*R - r*r - Δ*Δ) / (4 * r * Δ))
	if (l - l0).Sin() < 0 {
		C = -C
	}
	ω1 = (ω1 + C).Mod1()
	ω2 = (ω2 + C).Mod1()
	// Step 15.
	Δψ, Δε := nutation.Nutation(jde)
	ε := ε0 + Δε
	// Step 16.
	sε, cε := ε.Sincos()
	sα, cα := math.Sincos(α)
	α += .005693 * p * (cα*cl0*cε + sα*sl0) / cδ
	δ += .005693 * p * (cl0*cε*(sε/cε*cδ-sα*sδ) + cα*sδ*sl0)
	// Step 17.
	tδ := sδ / cδ
	Δα := (cε+sε*sα*tδ)*Δψ.Rad() - cα*tδ*Δε.Rad()
	Δδ := sε*cα*Δψ.Rad() + sα*Δε.Rad()
	αʹ := α + Δα
	δʹ := δ + Δδ
	sα0, cα0 := math.Sincos(α0)
	tδ0 := sδ0 / cδ0
	Δα0 := (cε+sε*sα0*tδ0)*Δψ.Rad() - cα0*tδ0*Δε.Rad()
	Δδ0 := sε*cα0*Δψ.Rad() + sα0*Δε.Rad()
	α0ʹ := α0 + Δα0
	δ0ʹ := δ0 + Δδ0
	// Step 18.
	sδʹ, cδʹ := math.Sincos(δʹ)
	sδ0ʹ, cδ0ʹ := math.Sincos(δ0ʹ)
	sα0ʹαʹ, cα0ʹαʹ := math.Sincos(α0ʹ - αʹ)
//.........这里部分代码省略.........

// Kepler4 returns an approximate solution to Kepler's equation.
//
// It is valid only for small values of e.
//
// Argument e is eccentricity, M is mean anomaly.
//
// Result E is eccentric anomaly.
func Kepler4(e float64, M unit.Angle) (E unit.Angle) {
	sm, cm := M.Sincos()
	return unit.Angle(math.Atan2(sm, cm-e)) // (30.8) p. 206
}

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

// RelativePosition returns the position angle of one body with respect to
// another.
//
// The position angle result is measured counter-clockwise from North.
func RelativePosition(r1, d1, r2, d2 unit.Angle) unit.Angle {
	sΔr, cΔr := (r2 - r1).Sincos()
	sd2, cd2 := d2.Sincos()
	return unit.Angle(math.Atan2(sΔr, cd2*d1.Tan()-sd2*cΔr))
}

// SunAltitude returns altitude of the Sun above the lunar horizon.
//
// Arguments η, θ are selenographic longitude and latitude of a site on the
// Moon, l0, b0 are selenographic coordinates of the Sun, as returned by
// Physical(), for example.
func SunAltitude(η, θ, l0, b0 unit.Angle) unit.Angle {
	c0 := math.Pi/2 - l0
	sb0, cb0 := b0.Sincos()
	sθ, cθ := θ.Sincos()
	return unit.Angle(math.Asin(sb0*sθ + cb0*cθ*(c0+η).Sin()))
}

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