第二十六章 太阳直角坐标 Rectangular Coordinates of the Sun

太阳的地心赤道直角坐标$X,Y,Z$。原点在地心。$X$轴的方向指向春分点（经度为0），$Y$轴平放在赤道面上， 经度是90°，$Z$轴方向是北极。

1. Date 平分点参考系

\begin{cases} X &= R\cos β\cos ☉\\[2ex] Y &= R(\cos β\sin ☉\cos ε - \sin β\sin ε)\\[2ex] Z &= R(\cos β\sin ☉\sin ε + \sin β\cos ε) \end{cases} 其中$☉,β,R$是太阳的地心几何黄经,黄纬,日地距离（可由solar.TrueVSOP87求出），$ε$是平黄赤交角。
因为 Date 黄道坐标中，太阳的纬度不超过 1.2 角秒， 所以可以看成$\cos β=1$

// Position returns rectangular coordinates referenced to the mean equinox
// of date.
// Date 平分点太阳地心直角坐标
func Position(e *pp.V87Planet, jde float64) (x, y, z float64) {
	// (26.1) p. 171
	s, β, R := solar.TrueVSOP87(e, jde)
	sε, cε := nutation.MeanObliquity(jde).Sincos()
	ss, cs := s.Sincos()
	sβ := β.Sin()
	x = R * cs
	y = R * (ss*cε - sβ*sε)
	z = R * (ss*sε + sβ*cε)
	return
}

2. J2000标准分点参考系和B1950平分点参考系

通过VSOP87计算给定时刻地球的相对于 J2000.0 分点的日心黄经$L$和黄纬$B$，及距离$R$ $$☉ = L + 180°， β = -B$$ 计算 \begin{cases} X &= R\cos β*cos ☉\\[2ex] Y &= R\cos β*sin ☉\\[2ex] Z &= R\sin β \end{cases}

上述坐标仍是黄道坐标，J2000参考系的赤道坐标为： \begin{cases} X_0 &= 1.000000000000X+0.000000440360Y-0.000000190919Z\\[2ex] Y_0 &= -0.000000479966X+0.917482137087Y-0.397776982902Z\\[2ex] Z_0 &= 0.000000000000X+0.397776982902Y+0.917482137087Z \end{cases}

// PositionJ2000 returns rectangular coordinates referenced to equinox J2000.
// J2000太阳直角坐标
func PositionJ2000(e *pp.V87Planet, jde float64) (x, y, z float64) {
	x, y, z = xyz(e, jde)
	// (26.3) p. 174
	return x + .00000044036*y - .000000190919*z,
		-.000000479966*x + .917482137087*y - .397776982902*z,
		.397776982902*y + .917482137087*z
}

func xyz(e *pp.V87Planet, jde float64) (x, y, z float64) {
	l, b, r := e.Position2000(jde)
	s := l + math.Pi
	β := -b
	ss, cs := s.Sincos()
	sβ, cβ := β.Sincos()
	// (26.2) p. 172
	x = r * cβ * cs
	y = r * cβ * ss
	z = r * sβ
	return
}

B1950平分点参考系坐标： \begin{cases} X_0 &=0.999925702634X+0.012189716217Y+0.000011134016Z\\[2ex] Y_0 &=-0.011179418036X+0.917413998946Y-0.397777041885Z\\[2ex] Z_0 &=-0.004859003787X+0.397747363646Y+0.917482111428Z \end{cases}

// PositionB1950 returns rectangular coordinates referenced to B1950.
// B1950 平分点参考系太阳直角坐标
//
// Results are referenced to the mean equator and equinox of the epoch B1950
// in the FK5 system, not FK4.
func PositionB1950(e *pp.V87Planet, jde float64) (x, y, z float64) {
	x, y, z = xyz(e, jde)
	return .999925702634*x + .012189716217*y + .000011134016*z,
		-.011179418036*x + .917413998946*y - .397777041885*z,
		-.004859003787*x + .397747363646*y + .917482111428*z
}

3. 任意其它平分点参考系

　　先计算出J2000标准分点标准$X_0,Y_0,Z_0$
　　然后按岁差计算中的方法计算出$ζ,z,θ$,
\begin{cases} X_x &= \cos ζ\cos z\cos θ - \sin ζ\sin z\\[2ex] X_y &= \sin ζ\cos z + \cos ζ\sin z\cos θ\\[2ex] X_z &= \cos ζ\sin θ\\[2ex] \end{cases} \begin{cases} Y_x &= -\cos ζ\sin z - \sin ζ\cos z\cos θ\\[2ex] Y_y &= \cos ζ\cos z - \sin ζ\sin z\cos θ\\[2ex] Y_z &= -\sin ζ\sin θ\\[2ex] \end{cases} \begin{cases} Z_x &= -\cos z\sin θ\\[2ex] Z_y &= -\sin z\sin θ\\[2ex] Z_z &= \cos θ \end{cases} 　　那么 \begin{cases} X′ &= X_xX_0 + Y_xY_0 + Z_xZ_0\\[2ex] Y′ &= X_yX_0 + Y_yY_0 + Z_yZ_0\\[2ex] Z′ &= X_zX_0 + Y_zY_0 + Z_zZ_0 \end{cases}

// PositionEquinox returns rectangular coordinates referenced to an arbitrary epoch.
// 任意其它平分点参考系太阳直角坐标
//
// Position will be computed for given Julian day "jde" but referenced to mean
// equinox "epoch" (year).
func PositionEquinox(e *pp.V87Planet, jde, epoch float64) (xp, yp, zp float64) {
	x0, y0, z0 := PositionJ2000(e, jde)
	t := (epoch - 2000) * .01
	ζ := base.Horner(t, ζt...) * t * math.Pi / 180 / 3600
	z := base.Horner(t, zt...) * t * math.Pi / 180 / 3600
	θ := base.Horner(t, θt...) * t * math.Pi / 180 / 3600
	sζ, cζ := math.Sincos(ζ)
	sz, cz := math.Sincos(z)
	sθ, cθ := math.Sincos(θ)
	xx := cζ*cz*cθ - sζ*sz
	xy := sζ*cz + cζ*sz*cθ
	xz := cζ * sθ
	yx := -cζ*sz - sζ*cz*cθ
	yy := cζ*cz - sζ*sz*cθ
	yz := -sζ * sθ
	zx := -cz * sθ
	zy := -sz * sθ
	zz := cθ
	return xx*x0 + yx*y0 + zx*z0,
		xy*x0 + yy*y0 + zy*z0,
		xz*x0 + yz*y0 + zz*z0
}