for Wolfgang

pre.go

package main

import (
    "fmt"
    "log"
    "math"

    "github.com/soniakeys/meeus/base"
    "github.com/soniakeys/meeus/coord"
    pp "github.com/soniakeys/meeus/planetposition"
    "github.com/soniakeys/meeus/precess"
    "github.com/soniakeys/meeus/solar"
    "github.com/soniakeys/sexagesimal"
    "github.com/soniakeys/unit"
)

func main() {
    show("Meeus:", coord.Ecliptic{
        unit.AngleFromDeg(46.865071),
        unit.AngleFromDeg(-2.539334),
    })
    show("Wolfgang:", coord.Ecliptic{
        unit.AngleFromDeg(46.865070344986776),
        unit.AngleFromDeg(-2.539333402296444),
    })
    show("Sonia:", ecl46a())
}

var jde = 2451439.50074

func show(who string, ecl coord.Ecliptic) {
    ep := base.JDEToJulianYear(jde)
    fmt.Printf("\n%-9s  %.2f  λ: %.12h  β: %.12h\n",
        who, ep, sexa.FmtAngle(ecl.Lon), sexa.FmtAngle(ecl.Lat))
    ep50 := base.JDEToJulianYear(base.B1950)
    precess.EclipticPosition(&ecl, &ecl, ep, ep50, 0, 0)
    fmt.Printf("Sonia:     %.2f  λ: %.12h  β: %.12h\n",
        ep50, sexa.FmtAngle(ecl.Lon), sexa.FmtAngle(ecl.Lat))
}

func ecl46a() coord.Ecliptic {
    earth, err := pp.LoadPlanet(pp.Earth)
    if err != nil {
        log.Fatal(err)
    }
    saturn, err := pp.LoadPlanet(pp.Saturn)
    if err != nil {
        log.Fatal(err)
    }
    s, β, R := solar.TrueVSOP87(earth, jde)
    ss, cs := s.Sincos()
    sβ := β.Sin()
    Δ := 9.
    var x, y, z float64
    var JDE float64
    f := func() {
        τ := base.LightTime(Δ)
        JDE = jde - τ
        l, b, r := saturn.Position(JDE)
        l, b = pp.ToFK5(l, b, JDE)
        sl, cl := l.Sincos()
        sb, cb := b.Sincos()
        x = r*cb*cl + R*cs
        y = r*cb*sl + R*ss
        z = r*sb + R*sβ
        Δ = math.Sqrt(x*x + y*y + z*z)
    }
    f()
    f()
    return coord.Ecliptic{
        unit.Angle(math.Atan2(y, x)),
        unit.Angle(math.Atan(z / math.Hypot(x, y))),
    }
}