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.

该提问来源于开源项目:soniakeys/meeus

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

11条回答 默认 最新

  • weixin_39861498 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 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 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
        return unit.Angle(math.Hypot((r2-r1).Rad()*cd1, (d2 - d1).Rad()))
    }
    
    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 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 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 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 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 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 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,
            445267.1114034, -.0018819, 1/545868, -1/113065000)).Mod1().Rad()
        M := unit.AngleFromDeg(base.Horner(T,
            357.5291092, 35999.0502909, -.0001535, 1/24490000)).Mod1().Rad()
        Mʹ := unit.AngleFromDeg(base.Horner(T, 134.9633964,
            477198.8675055, .0087414, 1/69699, -1/14712000)).Mod1().Rad()
        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 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 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)

    点赞 评论 复制链接分享

相关推荐