Golang math.Sincos函数代码示例

func newQs(JDE float64) *qs {
	var q qs
	q.t1 = JDE - 2411093
	q.t2 = q.t1 / 365.25
	q.t3 = (JDE-2433282.423)/365.25 + 1950
	q.t4 = JDE - 2411368
	q.t5 = q.t4 / 365.25
	q.t6 = JDE - 2415020
	q.t7 = q.t6 / 36525
	q.t8 = q.t6 / 365.25
	q.t9 = (JDE - 2442000.5) / 365.25
	q.t10 = JDE - 2409786
	q.t11 = q.t10 / 36525
	q.W0 = 5.095 * d * (q.t3 - 1866.39)
	q.W1 = 74.4*d + 32.39*d*q.t2
	q.W2 = 134.3*d + 92.62*d*q.t2
	q.W3 = 42*d - .5118*d*q.t5
	q.W4 = 276.59*d + .5118*d*q.t5
	q.W5 = 267.2635*d + 1222.1136*d*q.t7
	q.W6 = 175.4762*d + 1221.5515*d*q.t7
	q.W7 = 2.4891*d + .002435*d*q.t7
	q.W8 = 113.35*d - .2597*d*q.t7
	q.s1, q.c1 = math.Sincos(28.0817 * d)
	q.s2, q.c2 = math.Sincos(168.8112 * d)
	q.e1 = .05589 - .000346*q.t7
	q.sW0 = math.Sin(q.W0)
	q.s3W0 = math.Sin(3 * q.W0)
	q.s5W0 = math.Sin(5 * q.W0)
	q.sW1 = math.Sin(q.W1)
	q.sW2 = math.Sin(q.W2)
	q.sW3, q.cW3 = math.Sincos(q.W3)
	q.sW4, q.cW4 = math.Sincos(q.W4)
	q.sW7, q.cW7 = math.Sincos(q.W7)
	return &q
}

// AberrationRonVondrak uses the Ron-Vondrák expression to compute corrections
// due to aberration for equatorial coordinates of an object.
func AberrationRonVondrak(α, δ, jd float64) (Δα, Δδ float64) {
	T := base.J2000Century(jd)
	r := &rv{
		T:  T,
		L2: 3.1761467 + 1021.3285546*T,
		L3: 1.7534703 + 628.3075849*T,
		L4: 6.2034809 + 334.0612431*T,
		L5: 0.5995465 + 52.9690965*T,
		L6: 0.8740168 + 21.3299095*T,
		L7: 5.4812939 + 7.4781599*T,
		L8: 5.3118863 + 3.8133036*T,
		Lp: 3.8103444 + 8399.6847337*T,
		D:  5.1984667 + 7771.3771486*T,
		Mp: 2.3555559 + 8328.6914289*T,
		F:  1.6279052 + 8433.4661601*T,
	}
	var Xp, Yp, Zp float64
	// sum smaller terms first
	for i := 35; i >= 0; i-- {
		x, y, z := rvTerm[i](r)
		Xp += x
		Yp += y
		Zp += z
	}
	sα, cα := math.Sincos(α)
	sδ, cδ := math.Sincos(δ)
	// (23.4) p. 156
	return (Yp*cα - Xp*sα) / (c * cδ), -((Xp*cα+Yp*sα)*sδ - Zp*cδ) / c
}

// EqToEcl converts equatorial coordinates to ecliptic coordinates.
func (ecl *Ecliptic) EqToEcl(eq *Equatorial, ε *Obliquity) *Ecliptic {
	sα, cα := math.Sincos(eq.RA)
	sδ, cδ := math.Sincos(eq.Dec)
	ecl.Lon = math.Atan2(sα*ε.C+(sδ/cδ)*ε.S, cα) // (13.1) p. 93
	ecl.Lat = math.Asin(sδ*ε.C - cδ*ε.S*sα)      // (13.2) p. 93
	return ecl
}

// Horizontal computes data for a horizontal sundial.
//
// Argument φ is geographic latitude at which the sundial will be located,
// 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 Horizontal(φ, a float64) (lines []Line, center Point, u float64) {
	sφ, cφ := math.Sincos(φ)
	tφ := sφ / cφ
	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 := cφ*cH + sφ*tδ
			x := a * sH / Q
			y := a * (sφ*cH - cφ*tδ) / Q
			l.Points = append(l.Points, Point{x, y})
		}
		if len(l.Points) > 0 {
			lines = append(lines, l)
		}
	}
	center.Y = -a / tφ
	u = a / math.Abs(sφ)
	return
}

// 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, a float64) (lines []Line, center Point, u float64) {
	sφ, cφ := math.Sincos(φ)
	tφ := sφ / cφ
	sD, cD := math.Sincos(D)
	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
}

// ParallaxConstants computes parallax constants ρ sin φ′ and ρ cos φ′.
//
// Arguments are geographic latitude φ in radians and height h
// in meters above the ellipsoid.
func (e Ellipsoid) ParallaxConstants(φ, h float64) (s, c float64) {
	boa := 1 - e.Fl
	su, cu := math.Sincos(math.Atan(boa * math.Tan(φ)))
	s, c = math.Sincos(φ)
	hoa := h * 1e-3 / e.Er
	return su*boa + hoa*s, cu + hoa*c
}

// Andoyer computes approximate geodesic distance between two
// points p1 and p2 on the spheroid s.
// Computations use Andoyer-Lambert first order (in flattening) approximation.
//
// Reference: P.D. Thomas, Mathematical Models for Navigation Systems, TR-182,
// US Naval Oceanographic Office (1965).
func Andoyer(s *Spheroid, p1, p2 LL) float64 {
	a, f := s.A(), s.F()
	lat1, lon1 := p1.LatLon()
	lat2, lon2 := p2.LatLon()

	F, G := (lat1+lat2)/2.0, (lat1-lat2)/2.0
	L := math.Remainder(lon1-lon2, 360.0) / 2.0
	sinF, cosF := math.Sincos((math.Pi / 180.0) * F)
	sinG, cosG := math.Sincos((math.Pi / 180.0) * G)
	sinL, cosL := math.Sincos((math.Pi / 180.0) * L)

	S2 := sq(sinG*cosL) + sq(cosF*sinL)
	C2 := sq(cosG*cosL) + sq(sinF*sinL)
	S, C := math.Sqrt(S2), math.Sqrt(C2)

	omega := math.Atan2(S, C)
	R := (S * C) / omega
	D := 2.0 * omega * a
	H1 := (3.0*R - 1.0) / (2.0 * C2)
	H2 := (3.0*R + 1.0) / (2.0 * S2)
	dist := D * (1.0 + f*(H1*sq(sinF*cosG)-H2*sq(cosF*sinG)))

	if finite(dist) {
		return dist
	}

	// Antipodal points.
	if finite(R) {
		return 2.0 * s.Quad()
	}

	// Identical points.
	return 0.0
}

// TrueEquatorial returns the true geometric position of the Sun as equatorial coordinates.
func TrueEquatorial(jde float64) (α, δ float64) {
	s, _ := True(base.J2000Century(jde))
	ε := nutation.MeanObliquity(jde)
	ss, cs := math.Sincos(s)
	sε, cε := math.Sincos(ε)
	// (25.6, 25.7) p. 165
	return math.Atan2(cε*ss, cs), sε * ss
}

// 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 float64) float64 {
	sd2, cd2 := math.Sincos(d2)
	sr21, cr21 := math.Sincos(r2 - r1)
	sr32, cr32 := math.Sincos(r3 - r2)
	C1 := math.Atan2(sr21, cd2*math.Tan(d1)-sd2*cr21)
	C2 := math.Atan2(sr32, cd2*math.Tan(d3)-sd2*cr32)
	return C1 + C2
}

// 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 (eq *Equatorial) GalToEq(g *Galactic) *Equatorial {
	sdLon, cdLon := math.Sincos(g.Lon - galacticLon0)
	sgδ, cgδ := math.Sincos(galacticNorth.Dec)
	sb, cb := math.Sincos(g.Lat)
	y := math.Atan2(sdLon, cdLon*sgδ-(sb/cb)*cgδ)
	eq.RA = base.PMod(y+galacticNorth.RA, 2*math.Pi)
	eq.Dec = math.Asin(sb*sgδ + cb*cgδ*cdLon)
	return eq
}

// EqToGal converts equatorial coordinates to galactic coordinates.
//
// Equatorial coordinates must be referred to the standard equinox of B1950.0.
// For conversion to B1950, see package precess and utility functions in
// package "common".
func EqToGal(α, δ float64) (l, b float64) {
	sdα, cdα := math.Sincos(galacticNorth.RA - α)
	sgδ, cgδ := math.Sincos(galacticNorth.Dec)
	sδ, cδ := math.Sincos(δ)
	x := math.Atan2(sdα, cdα*sgδ-(sδ/cδ)*cgδ)       // (13.7) p. 94
	l = math.Mod(math.Pi+galacticLon0-x, 2*math.Pi) // (13.8) p. 94
	b = math.Asin(sδ*sgδ + cδ*cgδ*cdα)
	return
}

// EqToHz computes Horizontal coordinates from equatorial coordinates.
//
//	α: right ascension coordinate to transform, in radians
//	δ: declination coordinate to transform, in radians
//	φ: 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:
//
//	A: azimuth of observed point in radians, measured westward from the South.
//	h: elevation, or height of observed point in radians above horizon.
func EqToHz(α, δ, φ, ψ, st float64) (A, h float64) {
	H := sexa.Time(st).Rad() - ψ - α
	sH, cH := math.Sincos(H)
	sφ, cφ := math.Sincos(φ)
	sδ, cδ := math.Sincos(ψ)
	A = math.Atan2(sH, cH*sφ-(sδ/cδ)*cφ) // (13.5) p. 93
	h = math.Asin(sφ*sδ + cφ*cδ*cH)      // (13.6) p. 93
	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
}

// ApparentEquatorial returns the apparent position of the Sun as equatorial coordinates.
//
//	α: right ascension in radians
//	δ: declination in radians
func ApparentEquatorial(jde float64) (α, δ float64) {
	T := base.J2000Century(jde)
	λ := ApparentLongitude(T)
	ε := nutation.MeanObliquity(jde)
	sλ, cλ := math.Sincos(λ)
	// (25.8) p. 165
	sε, cε := math.Sincos(ε + .00256*math.Pi/180*math.Cos(node(T)))
	return math.Atan2(cε*sλ, cλ), math.Asin(sε * sλ)
}

// 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 float64) float64 {
	sd1, cd1 := math.Sincos(d1)
	sd2, cd2 := math.Sincos(d2)
	cd := sd1*sd2 + cd1*cd2*math.Cos(r1-r2) // (17.1) p. 109
	if cd < base.CosSmallAngle {
		return math.Acos(cd)
	}
	return math.Hypot((r2-r1)*cd1, d2-d1) // (17.2) p. 109
}