###### weixin_39861498
2020-12-02 13:00 阅读 1

# Seems typo in interp.go, julian.go, sidereal.go, angle.go, nutation.go, solstice.go, illum.go, moonposition.go, moonillum.go, moonphase.go, moon.go

Hi ,

First, thanks for your pretty library. This is really nice to go thru "Astronomical Algorithms" chapter by chapter with individual packages.

It seems that there's little mismatch between your code and the algorithm listed in the book when i was reading chapter 3 Interpolation with package `interp`.

algorithm in book (page 29): the third coeff is `3(H+J)`

code in interp.go:

``````go
func (d *Len5) Extremum() (x, y float64, err error) {
// (3.9) p. 29
nCoeff := []float64{
6*(d.b+d.c) - d.h - d.j,
0,
3 * (d.h + d.k),
2 * d.k,
}
den := d.k - 12*d.f
if den == 0 {
return 0, 0, ErrorExtremumOutside
}
n0, ok := iterate(0, func(n0 float64) float64 {
return base.Horner(n0, nCoeff...) / den
})
if !ok {
return 0, 0, ErrorNoConverge
}
if n0 < -2 || n0 > 2 {
return 0, 0, ErrorExtremumOutside
}
x = .5*d.xSum + .25*d.xDiff*n0
y = base.Horner(n0, d.interpCoeff...)
return x, y, nil
}
``````

the third coeff is `3 * (d.h + d.k)`, not `3 * (d.h + d.j)`

I guess this is a typo, pls kindly check.

• 点赞
• 写回答
• 关注问题
• 收藏
• 复制链接分享

#### 11条回答默认 最新

• weixin_39861498 2020-12-02 13:00

1 more typo found in `julian.go`, the JD of the beginning of Gregorian Calendar(1582-10-15 12:00:00) should be 2299161, while the current value is 2299151.

``````go
func JDToCalendar(jd float64) (year, month int, day float64) {
zf, f := math.Modf(jd + .5)
z := int64(zf)
a := z
if z >= 2299151 {
α := base.FloorDiv64(z*100-186721625, 3652425)
a = z + 1 + α - base.FloorDiv64(α, 4)
}
``````

I have created 1 pull request #16 for above changes, pls kindly review. Thanks

点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

1 more found in `sidereal.go`

``````go
var iau82 = []float64{24110.54841, 8640184.812866, 0.093104, 0.0000062}
``````

As per the IAU1982 Coeff, the 4th value should be -0.0000062, pls check. It looks like it's too small to impact the test cases listed in `sidereal_test.go`.

点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00
1. As per Formula 17.2(page 109), the δ is the average of the declinations of the two bodies. But in `angle.go`, the value is one of the declinations only.
``````go
func Sep(r1, d1, r2, d2 unit.Angle) unit.Angle {
sd1, cd1 := d1.Sincos()
sd2, cd2 := d2.Sincos()
cd := sd1*sd2 + cd1*cd2*(r1-r2).Cos() // (17.1) p. 109
if cd < base.CosSmallAngle {
return unit.Angle(math.Acos(cd))
}
// (17.2) p. 109
}
``````
1. As per the formula of Relative Position Angle(page 116), the Angle is measured counter-clockwise from the second body's North, thus Δr should be `r1-r2`. And I think it will impact the result if Δr is calculated as `r2-r1`.
``````go
func RelativePosition(r1, d1, r2, d2 unit.Angle) unit.Angle {
sΔr, cΔr := (r2 - r1).Sincos()
sd2, cd2 := d2.Sincos()
return unit.Angle(math.Atan2(sΔr, cd2*d1.Tan()-sd2*cΔr))
}
``````
点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

1 typo in `nutation.go`, as per formula in page 144, the 4th Coeff of `N` should be `1/56250`.

``````go
func Nutation(jde float64) (Δψ, Δε unit.Angle) {
T := base.J2000Century(jde)
D := base.Horner(T,
297.85036, 445267.11148, -0.0019142, 1./189474) * math.Pi / 180
M := base.Horner(T,
357.52772, 35999.050340, -0.0001603, -1./300000) * math.Pi / 180
N := base.Horner(T,
134.96298, 477198.867398, 0.0086972, 1./5620) * math.Pi / 180
``````
点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

the 3rd Coeff of June solstice calculation for years -1000 to 1000 is not correct: `solstice.go`

``````go
jc0 = []float64{1721233.25401, 365241.72562, -.05232, .00907, .00025}
``````
点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

in `illum.go`, the below formulas are not correct

``````go
func Venus84(r, Δ float64, i unit.Angle) float64 {
return base.Horner(i.Deg(), -4.4+5*math.Log10(r*Δ),
.0009, -.000239, .00000065)
}
...
func Saturn84(r, Δ float64, B, ΔU unit.Angle) float64 {
s := math.Abs(B.Sin())
return -8.88 + 5*math.Log10(r*Δ) + .044/math.Abs(ΔU.Deg()) -
2.6*s + 1.25*s*s
}
``````

should be

, fyi :)

点赞 评论 复制链接分享
• weixin_39890102 2020-12-02 13:00

~For `func Saturn()` the same applies as for `func Saturn84` in terms of considering `Math.abs(B.Sin())` and `(B.Sin())^2`.~

My bad. Forget this comment... `(B.Sin())^2 == Math.abs(B.Sin())^2`

点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

，should not be ```return -8.88 + 5*math.Log10(r*Δ) + .044*math.Abs(ΔU.Deg()) - 2.6*s + 1.25*s*s```?

点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

in `moonposition.go` & `moonillum.go`, the 3rd Coeff of sun mean anomaly calculation(page 338) is not correct.

``````go
func dmf(T float64) (D, M, Mʹ, F float64) {
D = base.Horner(T, 297.8501921*p, 445267.1114034*p,
-.0018819*p, p/545868, -p/113065000)
M = base.Horner(T, 357.5291092*p, 35999.0502909*p,
-.0001535*p, p/24490000)
Mʹ = base.Horner(T, 134.9633964*p, 477198.8675055*p,
.0087414*p, p/69699, -p/14712000)
F = base.Horner(T, 93.272095*p, 483202.0175233*p,
-.0036539*p, -p/3526000, p/863310000)
return
}
``````
``````go
func PhaseAngle3(jde float64) unit.Angle {
T := base.J2000Century(jde)
D := unit.AngleFromDeg(base.Horner(T, 297.8501921,
M := unit.AngleFromDeg(base.Horner(T,
Mʹ := unit.AngleFromDeg(base.Horner(T, 134.9633964,
return math.Pi - unit.Angle(D) + unit.AngleFromDeg(
-6.289*math.Sin(Mʹ)+
2.1*math.Sin(M)+
-1.274*math.Sin(2*D-Mʹ)+
-.658*math.Sin(2*D)+
-.214*math.Sin(2*Mʹ)+
-.11*math.Sin(D))
}
``````
点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

in `moonphase.go`, the 11th Coeff of correction for first/last quarter phase is not correct, should multiply one more `E`

``````go
// first or last corrections
func (m *mp) flc() float64 {
return -.62801*math.Sin(m.Mʹ) +
.17172*math.Sin(m.M)*m.E +
-.01183*math.Sin(m.Mʹ+m.M)*m.E +
.00862*math.Sin(2*m.Mʹ) +
.00804*math.Sin(2*m.F) +
.00454*math.Sin(m.Mʹ-m.M)*m.E +
.00204*math.Sin(2*m.M)*m.E*m.E +
-.0018*math.Sin(m.Mʹ-2*m.F) +
-.0007*math.Sin(m.Mʹ+2*m.F) +
-.0004*math.Sin(3*m.Mʹ) +
-.00034*math.Sin(2*m.Mʹ-m.M) +
.00032*math.Sin(m.M+2*m.F)*m.E +
.00032*math.Sin(m.M-2*m.F)*m.E +
-.00028*math.Sin(m.Mʹ+2*m.M)*m.E*m.E +
.00027*math.Sin(2*m.Mʹ+m.M)*m.E +
-.00017*math.Sin(m.Ω) +
-.00005*math.Sin(m.Mʹ-m.M-2*m.F) +
.00004*math.Sin(2*m.Mʹ+2*m.F) +
-.00004*math.Sin(m.Mʹ+m.M+2*m.F) +
.00004*math.Sin(m.Mʹ-2*m.M) +
.00003*math.Sin(m.Mʹ+m.M-2*m.F) +
.00003*math.Sin(3*m.M) +
.00002*math.Sin(2*m.Mʹ-2*m.F) +
.00002*math.Sin(m.Mʹ-m.M+2*m.F) +
-.00002*math.Sin(3*m.Mʹ+m.M)
}
``````
点赞 评论 复制链接分享
• weixin_39861498 2020-12-02 13:00

in `moon.go`, the following formula is not same as the one on page 376 to calculate the selenographic of sun.

``````go
func (m *moon) sun(λ, β unit.Angle, Δ float64, earth *pp.V87Planet) (l0, b0 unit.Angle) {
λ0, _, R := solar.ApparentVSOP87(earth, m.jde)
ΔR := unit.Angle(Δ / (R * base.AU))
λH := λ0 + math.Pi + ΔR.Mul(β.Cos()*(λ0-λ).Sin())
βH := ΔR * β
return m.lib(λH, βH)
}
``````

Also the 3rd Coeff of sun mean anomaly calculation is not correct(same issue as moonposition.go)

点赞 评论 复制链接分享