第二十五章 太阳坐标 Solar Coordinates

1. 低精度太阳黄经

当计算精度要求为0.01度，计算太阳位置时可假设地球运动是一个纯椭圆，也就说忽略月球及行星摄动。

• 太阳真黄经

  $T$为J2000起算的儒略世纪数
  \begin{align} 太阳几何平黄经：L_0 &= 280°.46646 + 36000°.76983T + 0°.0003032T^2\\[2ex] 太阳平近点角： M &= 357°.52911 + 35999°.05029T - 0°.0001537T^2\\[2ex] 地球轨道离心率： e &= 0.016708634 - 0.000042037T - 0.0000001267T^2\\[2ex] 太阳中心方程 ： C &= +(1°.914602 - 0°.004817T - 0°.000014T^2)\sin M\\[2ex] &+(0°.019993 - 0°.000101T)\sin 2M\\[2ex] &+ 0°.000289\sin 3M \end{align} 那么，太阳的真黄经是：$☉ = L_0 + C$ 真近点角是：$ν = M + C$
  日地距离(AU)为：$R =\frac {1.000001018(1-e^2)}{1+e\cos ν}$

  // True returns true geometric longitude and anomaly of the sun referenced to the mean equinox of date.
  // 计算太阳真黄经s，真近点角ν
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  //
  // Results:
  //  s = true geometric longitude, ☉
  //  ν = true anomaly
  func True(T float64) (s, ν unit.Angle) {
      // (25.2) p. 163
      L0 := unit.AngleFromDeg(base.Horner(T, 280.46646, 36000.76983, 0.0003032))
      M := MeanAnomaly(T)
      C := unit.AngleFromDeg(base.Horner(T, 1.914602, -0.004817, -.000014)*
          M.Sin() +
          (0.019993-.000101*T)*M.Mul(2).Sin() +
          0.000289*M.Mul(3).Sin())
      return (L0 + C).Mod1(), (M + C).Mod1()
  }
  
  // MeanAnomaly returns the mean anomaly of Earth at the given T.
  // 太阳平近点角
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  //
  // Result is not normalized to the range 0..2π.
  func MeanAnomaly(T float64) unit.Angle {
      // (25.3) p. 163
      return unit.AngleFromDeg(base.Horner(T, 357.52911, 35999.05029, -0.0001537))
  }
  
  // Eccentricity returns eccentricity of the Earth's orbit around the sun.
  // 地球轨道离心率
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  func Eccentricity(T float64) float64 {
      // (25.4) p. 163
      return base.Horner(T, 0.016708634, -0.000042037, -0.0000001267)
  }
  
  // Radius returns the Sun-Earth distance in AU.
  // 日地距离，单位为 AU
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  func Radius(T float64) float64 {
      _, ν := True(T)
      e := Eccentricity(T)
      // (25.5) p. 164
      return 1.000001018 * (1 - e*e) / (1 + e*ν.Cos())
  }
  

• 太阳视黄经

  太阳视黄经$λ$ = 太阳真黄经$☉$ + 章动修正 + 光行差修正

  如果精度要求不高，可以采用以下公式： \begin{cases} Ω &= 125°.04 - 1934°.136T\\[2ex] λ &= ☉ - 0°.00569 -0°.00478\sin Ω \end{cases}

  // ApparentLongitude returns apparent longitude of the Sun referenced
  // to the true equinox of date.
  // 太阳视黄经，考虑了章动和光行差
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  //
  // Result includes correction for nutation and aberration.
  func ApparentLongitude(T float64) unit.Angle {
      Ω := node(T)
      s, _ := True(T)
      return s - unit.AngleFromDeg(.00569) -
          unit.AngleFromDeg(.00478).Mul(Ω.Sin())
  }
  
  func node(T float64) unit.Angle {
      return unit.AngleFromDeg(125.04 - 1934.136*T)
  }
  

  J2000的太阳真黄经$☉_2000 = ☉ - 0°.01397(year-2000), 1900\leq year \leq 2100$

  // True2000 returns true geometric longitude and anomaly of the sun referenced to equinox J2000.
  //
  // Argument T is the number of Julian centuries since J2000.
  // See base.J2000Century.
  //
  // Results are accurate to .01 degree for years 1900 to 2100.
  //
  // Results:
  //  s = true geometric longitude, ☉
  //  ν = true anomaly
  // J2000的太阳真黄经，真近点角
  func True2000(T float64) (s, ν unit.Angle) {
      s, ν = True(T)
      s -= unit.AngleFromDeg(.01397).Mul(T * 100)
      return
  }
  

• 太阳地心赤经$α$，地心赤纬$δ$

  \begin{cases} \tan α = \frac {\cos ε\sin ☉}{\cos ☉}\\[2ex] \sin δ = \sin ε\sin ☉ \end{cases}

  // TrueEquatorial returns the true geometric position of the Sun as equatorial coordinates.
  // 太阳真赤经，真赤纬
  func TrueEquatorial(jde float64) (α unit.RA, δ unit.Angle) {
      s, _ := True(base.J2000Century(jde))
      ε := nutation.MeanObliquity(jde)
      ss, cs := s.Sincos()
      sε, cε := ε.Sincos()
      // (25.6, 25.7) p. 165
      α = unit.RAFromRad(math.Atan2(cε*ss, cs))
      δ = unit.Angle(math.Asin(sε * ss))
      return
  }
  

• 太阳地心视赤经，视赤纬

  $☉$补上黄经章动及光行差得到太阳视黄经$λ$，$ε$补上交角章动$+0.00256\cos Ω$ \begin{cases} \tan α = \frac {\cos ε\sin λ}{\cos λ}\\[2ex] \sin δ = \sin ε\sin λ \end{cases}

  // ApparentEquatorial returns the apparent position of the Sun as equatorial coordinates.
  // 太阳视赤经，视赤纬
  //
  //  α: right ascension in radians
  //  δ: declination in radians
  func ApparentEquatorial(jde float64) (α unit.RA, δ unit.Angle) {
      T := base.J2000Century(jde)
      λ := ApparentLongitude(T)
      ε := nutation.MeanObliquity(jde)
      sλ, cλ := λ.Sincos()
      // (25.8) p. 165
      ε += unit.AngleFromDeg(.00256).Mul(node(T).Cos())
      sε, cε := ε.Sincos()
      α = unit.RAFromRad(math.Atan2(cε*sλ, cλ))
      δ = unit.Angle(math.Asin(sε * sλ))
      return
  }
  

2. 高精度太阳黄经

• 太阳地心黄经，地心黄纬 先采用VSOP87理论计算地球位置得到地球的日心黄经$L$，黄纬$β$，日地距离$R$ 再计算地心黄经$☉ = L + 180°$,黄纬$β=-B$
  令$λ′ = ☉ - 1°.397T - 0°.00031T^2$
  那么 \begin{cases} Δ☉ = -0”.09033\\[2ex] Δβ = +0”.03916(\cos(λ′) - \sin(λ′)) \end{cases}

  // TrueVSOP87 returns the true geometric position of the sun as ecliptic coordinates.
  // 根据VSOP87理论计算太阳真黄经，真黄纬，日地距离
  //
  // Result computed by full VSOP87 theory.  Result is at equator and equinox
  // of date in the FK5 frame.  It does not include nutation or aberration.
  //
  //  s: ecliptic longitude
  //  β: ecliptic latitude
  //  R: range in AU
  func TrueVSOP87(e *pp.V87Planet, jde float64) (s, β unit.Angle, R float64) {
      l, b, r := e.Position(jde) //VSOP87算出的地球的日心黄经，黄纬，日地距离
      s = l + math.Pi
      // FK5 correction.
      λp := base.Horner(base.J2000Century(jde),
          s.Rad(), -1.397*math.Pi/180, -.00031*math.Pi/180)
      sλp, cλp := math.Sincos(λp)
      Δβ := unit.AngleFromSec(.03916).Mul(cλp - sλp)
      // (25.9) p. 166
      s -= unit.AngleFromSec(.09033)
      return s.Mod1(), Δβ - b, r
  }
  

• 太阳地心视黄经

  与日心视黄经一样，对日心黄经进行章动和光行差修正 章动可以由章动公式求得，光行差可以由以下公式求得$$-20”.4898/R$$

  // ApparentVSOP87 returns the apparent position of the sun as ecliptic coordinates.
  // 根据VSOP87理论计算太阳视黄经，视黄纬，日地距离，考虑了章动和光行差
  // 即真黄经+黄经章动+光行差，真黄纬，日地距离不变
  //
  // Result computed by VSOP87, at equator and equinox of date in the FK5 frame,
  // and includes effects of nutation and aberration.
  //
  //  λ: ecliptic longitude
  //  β: ecliptic latitude
  //  R: range in AU
  func ApparentVSOP87(e *pp.V87Planet, jde float64) (λ, β unit.Angle, R float64) {
      // note: see duplicated code in ApparentEquatorialVSOP87.
      s, β, R := TrueVSOP87(e, jde)
      Δψ, _ := nutation.Nutation(jde)
      a := aberration(R)
      return s + Δψ + a, β, R
  }
  
  // Low precision formula.  The high precision formula is not implemented
  // because the low precision formula already gives position results to the
  // accuracy given on p. 165.  The high precision formula the represents lots
  // of typing with associated chance of typos, and no way to test the result.
  // 低精度光行差修正项
  func aberration(R float64) unit.Angle {
      // (25.10) p. 167
      return unit.AngleFromSec(-20.4898).Div(R)
  }
  

• 太阳地心视赤经，视赤纬

  算法同低精度太阳地心视赤经，视赤纬一致

  // ApparentEquatorialVSOP87 returns the apparent position of the sun as equatorial coordinates.
  // 根据VSOP87理论计算太阳视赤经，视赤纬，日地距离，考虑了章动和光行差
  // 即先计算视黄经，视黄纬，此时考虑交角章动，用真黄赤交角转为赤道坐标
  //
  // Result computed by VSOP87, at equator and equinox of date in the FK5 frame,
  // and includes effects of nutation and aberration.
  //
  //  α: right ascension
  //  δ: declination
  //  R: range in AU
  func ApparentEquatorialVSOP87(e *pp.V87Planet, jde float64) (α unit.RA, δ unit.Angle, R float64) {
      // note: duplicate code from ApparentVSOP87 so we can keep Δε.
      // see also duplicate code in time.E().
      s, β, R := TrueVSOP87(e, jde)
      Δψ, Δε := nutation.Nutation(jde)
      a := aberration(R)
      λ := s + Δψ + a
      ε := nutation.MeanObliquity(jde) + Δε
      sε, cε := ε.Sincos()
      α, δ = coord.EclToEq(λ, β, sε, cε)
      return
  }