第二十九章 Ephemeris for Physical Observations of the Sun

1. 日面计算

$P$ 为日轴方位角，自日面北点量起的太阳自转轴北端的方位角，向东为正
$B_0$ 为日面中心的日面纬度
$L_0$ 为日面中心的日面经度

根据之前章节的算法，先计算太阳视黄经$λ$(包含光行差修正)，$λ’$为$λ$的黄经章动修正值，$ε$为真黄赤交角
再计算 \begin{cases} θ = (JD - 2398220)\frac {360°}{25.38}\\[2ex] I = 7°.25,(太阳赤道与黄道的倾角)\\[2ex] K = 73°.6667 + 1°.3958333\frac {JD - 2396758}{36525},(太阳赤道与黄道的升交点的黄经)\\[2ex] \tan x = -\cos λ’\tan ε\\[2ex] \tan y = -\cos (λ - K)\tan I \end{cases} 那么 \begin{cases} P = x+y\\[2ex] \sin B_0 = \sin(λ - K)\sin I\\[2ex] \tan η = \frac {-\sin(λ - K)\cos I}{-\cos (λ - K)}\\[2ex] L_0 = η-θ \end{cases}

// Ephemeris returns the apparent orientation of the sun at the given jd.
// 计算日面参数
//
// Results:
//	P:  Position angle of the solar north pole.
//	B0: Heliographic latitude of the center of the solar disk.
//	L0: Heliographic longitude of the center of the solar disk.
func Ephemeris(jd float64, e *pp.V87Planet) (P, B0, L0 unit.Angle) {
	θ := unit.Angle((jd - 2398220) * 2 * math.Pi / 25.38)
	I := unit.AngleFromDeg(7.25)
	K := unit.AngleFromDeg(73.6667) +
		unit.AngleFromDeg(1.3958333).Mul((jd-2396758)/base.JulianCentury)

	L, _, R := solar.TrueVSOP87(e, jd)
	Δψ, Δε := nutation.Nutation(jd)
	ε0 := nutation.MeanObliquity(jd)
	ε := ε0 + Δε
	λ := L - unit.AngleFromSec(20.4898).Div(R)
	λp := λ + Δψ

	sλK, cλK := (λ - K).Sincos()
	sI, cI := I.Sincos()

	tx := -(λp.Cos() * ε.Tan())
	ty := -(cλK * I.Tan())
	P = unit.Angle(math.Atan(tx) + math.Atan(ty))
	B0 = unit.Angle(math.Asin(sλK * sI))
	η := unit.Angle(math.Atan2(-sλK*cI, -cλK))
	L0 = (η - θ).Mod1()
	return
}

// Cycle returns the jd of the start of the given synodic rotation.
//
// Argument c is the "Carrington" cycle number.
//
// Result is a dynamical time (not UT).
func Cycle(c int) (jde float64) {
	cf := float64(c)
	jde = 2398140.227 + 27.2752316*cf
	m := 281.96*math.Pi/180 + 26.882476*math.Pi/180*cf
	s2m, c2m := math.Sincos(2 * m)
	return jde + .1454*math.Sin(m) - .0085*s2m - .0141*c2m
}