Golang base.JDEToJulianYear函数代码示例

// Exercise, p. 136.
func TestPosition(t *testing.T) {
	eqFrom := &coord.Equatorial{
		base.NewRA(2, 31, 48.704).Rad(),
		base.NewAngle(false, 89, 15, 50.72).Rad(),
	}
	eqTo := &coord.Equatorial{}
	mα := base.NewHourAngle(false, 0, 0, 0.19877)
	mδ := base.NewAngle(true, 0, 0, 0.0152)
	for _, tc := range []struct {
		α, δ string
		jde  float64
	}{
		{"1 22 33.90", "88 46 26.18", base.BesselianYearToJDE(1900)},
		{"3 48 16.43", "89 27 15.38", base.JulianYearToJDE(2050)},
		{"5 53 29.17", "89 32 22.18", base.JulianYearToJDE(2100)},
	} {
		epochTo := base.JDEToJulianYear(tc.jde)
		precess.Position(eqFrom, eqTo, 2000.0, epochTo, mα, mδ)
		αStr := fmt.Sprintf("%.2x", base.NewFmtRA(eqTo.RA))
		δStr := fmt.Sprintf("%.2x", base.NewFmtAngle(eqTo.Dec))
		if αStr != tc.α {
			t.Fatal("got:", αStr, "want:", tc.α)
		}
		if δStr != tc.δ {
			t.Fatal(δStr)
		}
	}
}

func TestPrecessor_Precess(t *testing.T) {
	// Exercise, p. 136.
	eqFrom := &coord.Equatorial{
		RA:  base.NewRA(2, 31, 48.704).Rad(),
		Dec: base.NewAngle(false, 89, 15, 50.72).Rad(),
	}
	mα := base.NewHourAngle(false, 0, 0, .19877)
	mδ := base.NewAngle(false, 0, 0, -.0152)
	epochs := []float64{
		base.JDEToJulianYear(base.B1900),
		2050,
		2100,
	}
	answer := []string{
		"α = 1ʰ22ᵐ33ˢ.90   δ = +88°46′26″.18",
		"α = 3ʰ48ᵐ16ˢ.43   δ = +89°27′15″.38",
		"α = 5ʰ53ᵐ29ˢ.17   δ = +89°32′22″.18",
	}
	eqTo := &coord.Equatorial{}
	for i, epochTo := range epochs {
		precess.Position(eqFrom, eqTo, 2000, epochTo, mα, mδ)
		if answer[i] != fmt.Sprintf("α = %0.2d   δ = %+0.2d",
			base.NewFmtRA(eqTo.RA), base.NewFmtAngle(eqTo.Dec)) {
			t.Fatal(i)
		}
	}
}

// Exercise, p. 136.
func TestPosition(t *testing.T) {
	eqFrom := &coord.Equatorial{
		unit.NewRA(2, 31, 48.704),
		unit.NewAngle(' ', 89, 15, 50.72),
	}
	eqTo := &coord.Equatorial{}
	mα := unit.HourAngleFromSec(0.19877)
	mδ := unit.AngleFromSec(-0.0152)
	for _, tc := range []struct {
		α, δ string
		jde  float64
	}{
		{"1ʰ22ᵐ33.90ˢ", "88°46′26.18″", base.BesselianYearToJDE(1900)},
		{"3ʰ48ᵐ16.43ˢ", "89°27′15.38″", base.JulianYearToJDE(2050)},
		{"5ʰ53ᵐ29.17ˢ", "89°32′22.18″", base.JulianYearToJDE(2100)},
	} {
		epochTo := base.JDEToJulianYear(tc.jde)
		precess.Position(eqFrom, eqTo, 2000.0, epochTo, mα, mδ)
		αStr := fmt.Sprintf("%.2s", sexa.FmtRA(eqTo.RA))
		δStr := fmt.Sprintf("%.2s", sexa.FmtAngle(eqTo.Dec))
		if αStr != tc.α {
			t.Fatal("got:", αStr, "want:", tc.α)
		}
		if δStr != tc.δ {
			t.Fatal(δStr)
		}
	}
}

func Test255(t *testing.T) {
	for _, d := range td {
		_, f := math.Modf(.5 + d.f(base.JDEToJulianYear(d.jNom)))
		if int(f*24+.5) != d.hour {
			t.Errorf("got %d, expected %d", int(f*24+.5), d.hour)
		}
	}
}

func ExampleEclipticPrecessor_ReduceElements() {
	// Example 24.a, p. 160.
	ele := &elementequinox.Elements{
		Inc:  47.122 * math.Pi / 180,
		Peri: 151.4486 * math.Pi / 180,
		Node: 45.7481 * math.Pi / 180,
	}
	JFrom := base.JDEToJulianYear(base.BesselianYearToJDE(1744))
	JTo := base.JDEToJulianYear(base.BesselianYearToJDE(1950))
	p := precess.NewEclipticPrecessor(JFrom, JTo)
	p.ReduceElements(ele, ele)
	fmt.Printf("i = %.4f\n", ele.Inc*180/math.Pi)
	fmt.Printf("Ω = %.4f\n", ele.Node*180/math.Pi)
	fmt.Printf("ω = %.4f\n", ele.Peri*180/math.Pi)
	// Output:
	// i = 47.1380
	// Ω = 48.6037
	// ω = 151.4782
}

func ExamplePoly1900to1997() {
	// Example 10.a, p. 78.
	jd := julian.TimeToJD(time.Date(1977, 2, 18, 3, 37, 40, 0, time.UTC))
	year := base.JDEToJulianYear(jd)
	fmt.Printf("julian year %.1f\n", year)
	fmt.Printf("%+.1f seconds\n", deltat.Poly1900to1997(jd))
	// Output:
	// julian year 1977.1
	// +47.1 seconds
}

func ExampleEclipticPrecessor_ReduceElements() {
	// Example 24.a, p. 160.
	ele := &elementequinox.Elements{
		Inc:  unit.AngleFromDeg(47.122),
		Peri: unit.AngleFromDeg(151.4486),
		Node: unit.AngleFromDeg(45.7481),
	}
	JFrom := base.JDEToJulianYear(base.BesselianYearToJDE(1744))
	JTo := base.JDEToJulianYear(base.BesselianYearToJDE(1950))
	p := precess.NewEclipticPrecessor(JFrom, JTo)
	p.ReduceElements(ele, ele)
	fmt.Printf("i = %.4f\n", ele.Inc.Deg())
	fmt.Printf("Ω = %.4f\n", ele.Node.Deg())
	fmt.Printf("ω = %.4f\n", ele.Peri.Deg())
	// Output:
	// i = 47.1380
	// Ω = 48.6037
	// ω = 151.4782
}

// Position returns ecliptic position of planets at equinox and ecliptic of date.
//
// Argument jde is the date for which positions are desired.
//
// Results are positions consistent with those from Meeus's Apendix III,
// that is, at equinox and ecliptic of date.
//
//  L is heliocentric longitude in radians.
//  B is heliocentric latitude in radians.
//  R is heliocentric range in AU.
func (vt *V87Planet) Position(jde float64) (L, B, R float64) {
	L, B, R = vt.Position2000(jde)
	eclFrom := &coord.Ecliptic{
		Lat: B,
		Lon: L,
	}
	eclTo := &coord.Ecliptic{}
	epochFrom := 2000.0
	epochTo := base.JDEToJulianYear(jde)
	precess.EclipticPosition(eclFrom, eclTo, epochFrom, epochTo, 0, 0)
	return eclTo.Lon, eclTo.Lat, R
}

func ExampleEclipticPosition() {
	// Example 21.c, p. 137.
	eclFrom := &coord.Ecliptic{
		Lat: 1.76549 * math.Pi / 180,
		Lon: 149.48194 * math.Pi / 180,
	}
	eclTo := &coord.Ecliptic{}
	epochFrom := 2000.0
	epochTo := base.JDEToJulianYear(julian.CalendarJulianToJD(-214, 6, 30))
	precess.EclipticPosition(eclFrom, eclTo, epochFrom, epochTo, 0, 0)
	fmt.Printf("%.3f\n", eclTo.Lon*180/math.Pi)
	fmt.Printf("%+.3f\n", eclTo.Lat*180/math.Pi)
	// Output:
	// 118.704
	// +1.615
}

func ExampleEclipticPosition() {
	// Example 21.c, p. 137.
	eclFrom := &coord.Ecliptic{
		Lat: unit.AngleFromDeg(1.76549),
		Lon: unit.AngleFromDeg(149.48194),
	}
	eclTo := &coord.Ecliptic{}
	epochFrom := 2000.0
	epochTo := base.JDEToJulianYear(julian.CalendarJulianToJD(-214, 6, 30))
	precess.EclipticPosition(eclFrom, eclTo, epochFrom, epochTo, 0, 0)
	fmt.Printf("%.3f\n", eclTo.Lon.Deg())
	fmt.Printf("%+.3f\n", eclTo.Lat.Deg())
	// Output:
	// 118.704
	// +1.615
}

func ExamplePositionRonVondrak() {
	// Example 23.b, p. 156
	jd := julian.CalendarGregorianToJD(2028, 11, 13.19)
	eq := &coord.Equatorial{
		RA:  unit.NewRA(2, 44, 11.986),
		Dec: unit.NewAngle(' ', 49, 13, 42.48),
	}
	apparent.PositionRonVondrak(eq, eq, base.JDEToJulianYear(jd),
		unit.HourAngleFromSec(.03425),
		unit.AngleFromSec(-.0895))
	fmt.Printf("α = %0.3d\n", sexa.FmtRA(eq.RA))
	fmt.Printf("δ = %0.2d\n", sexa.FmtAngle(eq.Dec))
	// Output:
	// α = 2ʰ46ᵐ14ˢ.392
	// δ = 49°21′07″.45
}

func ExamplePosition() {
	// Example 23.a, p. 152
	jd := julian.CalendarGregorianToJD(2028, 11, 13.19)
	eq := &coord.Equatorial{
		sexa.NewRA(2, 44, 11.986).Rad(),
		sexa.NewAngle(false, 49, 13, 42.48).Rad(),
	}
	apparent.Position(eq, eq, 2000, base.JDEToJulianYear(jd),
		sexa.NewHourAngle(false, 0, 0, 0.03425),
		sexa.NewAngle(true, 0, 0, 0.0895))
	fmt.Printf("α = %0.3d\n", sexa.NewFmtRA(eq.RA))
	fmt.Printf("δ = %0.2d\n", sexa.NewFmtAngle(eq.Dec))
	// Output:
	// α = 2ʰ46ᵐ14ˢ.390
	// δ = 49°21′07″.45
}

func ExamplePosition() {
	// Example 21.b, p. 135.
	eq := &coord.Equatorial{
		base.NewRA(2, 44, 11.986).Rad(),
		base.NewAngle(false, 49, 13, 42.48).Rad(),
	}
	epochFrom := 2000.0
	jdTo := julian.CalendarGregorianToJD(2028, 11, 13.19)
	epochTo := base.JDEToJulianYear(jdTo)
	precess.Position(eq, eq, epochFrom, epochTo,
		base.NewHourAngle(false, 0, 0, 0.03425),
		base.NewAngle(true, 0, 0, 0.0895))
	fmt.Printf("%0.3d\n", base.NewFmtRA(eq.RA))
	fmt.Printf("%+0.2d\n", base.NewFmtAngle(eq.Dec))
	// Output:
	// 2ʰ46ᵐ11ˢ.331
	// +49°20′54″.54
}

// Positions returns positions of the eight major moons of Saturn.
//
// Results returned in argument pos, which must not be nil.
//
// Result units are Saturn radii.
func Positions(jde float64, earth, saturn *pp.V87Planet, pos *[8]XY) {
	s, β, R := solar.TrueVSOP87(earth, jde)
	ss, cs := s.Sincos()
	sβ := β.Sin()
	Δ := 9.
	var x, y, z float64
	var JDE float64
	f := func() {
		τ := base.LightTime(Δ)
		JDE = jde - τ
		l, b, r := saturn.Position(JDE)
		l, b = pp.ToFK5(l, b, JDE)
		sl, cl := l.Sincos()
		sb, cb := b.Sincos()
		x = r*cb*cl + R*cs
		y = r*cb*sl + R*ss
		z = r*sb + R*sβ
		Δ = math.Sqrt(x*x + y*y + z*z)
	}
	f()
	f()
	λ0 := unit.Angle(math.Atan2(y, x))
	β0 := unit.Angle(math.Atan(z / math.Hypot(x, y)))
	ecl := &coord.Ecliptic{λ0, β0}
	precess.EclipticPosition(ecl, ecl,
		base.JDEToJulianYear(jde), base.JDEToJulianYear(base.B1950), 0, 0)
	λ0, β0 = ecl.Lon, ecl.Lat
	q := newQs(JDE)
	s4 := [9]r4{{}, // 0 unused
		q.mimas(),
		q.enceladus(),
		q.tethys(),
		q.dione(),
		q.rhea(),
		q.titan(),
		q.hyperion(),
		q.iapetus(),
	}
	var X, Y, Z [9]float64
	for j := 1; j <= 8; j++ {
		u := s4[j].λ - s4[j].Ω
		w := s4[j].Ω - 168.8112*d
		su, cu := math.Sincos(u)
		sw, cw := math.Sincos(w)
		sγ, cγ := math.Sincos(s4[j].γ)
		r := s4[j].r
		X[j] = r * (cu*cw - su*cγ*sw)
		Y[j] = r * (su*cw*cγ + cu*sw)
		Z[j] = r * su * sγ
	}
	Z[0] = 1
	sλ0, cλ0 := λ0.Sincos()
	sβ0, cβ0 := β0.Sincos()
	var A, B, C [9]float64
	for j := range X {
		a := X[j]
		b := q.c1*Y[j] - q.s1*Z[j]
		c := q.s1*Y[j] + q.c1*Z[j]
		a, b =
			q.c2*a-q.s2*b,
			q.s2*a+q.c2*b
		A[j], b =
			a*sλ0-b*cλ0,
			a*cλ0+b*sλ0
		B[j], C[j] =
			b*cβ0+c*sβ0,
			c*cβ0-b*sβ0
	}
	D := math.Atan2(A[0], C[0])
	sD, cD := math.Sincos(D)
	for j := 1; j <= 8; j++ {
		X[j] = A[j]*cD - C[j]*sD
		Y[j] = A[j]*sD + C[j]*cD
		Z[j] = B[j]
		d := X[j] / s4[j].r
		X[j] += math.Abs(Z[j]) / k[j] * math.Sqrt(1-d*d)
		W := Δ / (Δ + Z[j]/2475)
		pos[j-1].X = X[j] * W
		pos[j-1].Y = Y[j] * W
	}
	return
}