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

[mooncaker816]

Hi @soniakeys ,

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:

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.

[mooncaker816]

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.

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

[mooncaker816]

1 more found in sidereal.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.

[mooncaker816]

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.

   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()))
   }

2. 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.

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))
}

[mooncaker816]

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

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

[mooncaker816]

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

jc0 = []float64{1721233.25401, 365241.72562, -.05232, .00907, .00025}

[mooncaker816]

in illum.go, the below formulas are not correct

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

@commenthol , fyi :)

[commenthol]

~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

[mooncaker816]

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

[mooncaker816]

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

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
}

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))
}

[mooncaker816]

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

// 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)
}

[mooncaker816]

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

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)