How do I use ECIToLookAngles correctly?

[ThomasHabets]

I made a test program to give me long and lat correctly. They seem to line up with what GPredict shows me. But I fail to convert these to LookAngles.

The test program below should say that the elevation is straight up, but instead the output of this code is:

ISS pos: 51.75 0.15
LookAngles: 1.32 -0.05 2631.53

Any ideas? Is this using jday incorrectly?

package main

import (
        "flag"
        "fmt"
        "log"
        "math"
        "time"

        "github.com/joshuaferrara/go-satellite"
)

func ymd(ts time.Time) (int, int, int, int, int, int) {
        y, mm, d := ts.Date()
        h, m, s := ts.Hour(), ts.Minute(), ts.Second()
        return y, int(mm), d, h, m, s
}

func gstimeFromDate(ts time.Time) float64 {
        return satellite.GSTimeFromDate(ymd(ts))
}

func fixLong(in float64) float64 {
        for in < 0 {
                in += 2 * math.Pi
        }
        for in > 2*math.Pi {
                in -= 2 * math.Pi
        }
        return in * 180 / math.Pi
}

func main() {
        flag.Parse()

        // Set up ISS TLE.                                                                                                                                            
tle1 := "1 25544U 98067A   21290.96791059  .00007152  00000-0  13913-3 0  9995"
        tle2 := "2 25544  51.6432  95.6210 0004029 117.2302  27.3327 15.48732973307628"
        sat := satellite.TLEToSat(tle1, tle2, "wgs84")

        // Pick a time when ISS is over London, UK.                                                                                                                   ts, err := time.Parse("2006-01-02 15:04:05", "2021-10-20 09:51:06")
        if err != nil {
                log.Fatal(err)
        }
        y, mm, d, h, m, s := ymd(ts)

        // Find ISS pos.                                                                                                                                              pos, _ := satellite.Propagate(sat, y, mm, d, h, m, s)

        // Find LLA for ISS.                                                                                                                                          gmst := satellite.GSTimeFromDate(ymd(ts))
        _, _, latlong := satellite.ECIToLLA(pos, gmst)
        fmt.Printf("ISS pos: %.2f %.2f\n",
                latlong.Latitude*180.0/math.Pi,
                fixLong(latlong.Longitude))

        // Now calculate LookAngles for *exactly* underneath the ISS.                                                                                                 jday := satellite.JDay(y, mm, d, h, m, s)
        ang := satellite.ECIToLookAngles(pos, satellite.LatLong{
                Latitude:  latlong.Latitude,
                Longitude: latlong.Longitude,
        }, 0, jday-2451545.0)
        fmt.Printf("LookAngles: %.2f %.2f %.2f\n", ang.Az, ang.El, ang.Rg)
}

[bclswl0827]

Hello, here is a correct example :)

package main

import (
    "fmt"
    "math"
    "time"

    "github.com/joshuaferrara/go-satellite"
)

func main() {
    // declare TLE content
    tle1 := "1 43823U 18100A   22156.97857081 -.00000333  00000+0  00000+0 0  9999"
    tle2 := "2 43823   0.0538 129.7336 0002248 306.5338 298.6476  1.00271874 12907"
    // declare current time
    year, month, day, hour, minute, second := time.Now().Year(),
        int(time.Now().Month()),
        time.Now().Day(),
        time.Now().Hour(),
        time.Now().Minute(),
        time.Now().Second()
    // initialize satellite
    sat := satellite.TLEToSat(tle1, tle2, "wgs84")
    // get the satellite position
    position, _ := satellite.Propagate(
        sat, year, month, day,
        hour, minute, second,
    )
    // convert the current time to Galileo system time (GST)
    gst := satellite.GSTimeFromDate(
        year, month, day,
        hour, minute, second,
    )
    // get satellite coordinates in radians
    _, _, latlng := satellite.ECIToLLA(position, gst)
    // declare my current location, altitude
    location, altitude := satellite.LatLong{}, 0.0
    // all latitude, longitude and azimuth angles need to be converted to radians
    location.Latitude = 29.921561 * math.Pi / 180
    location.Longitude = 106.629555 * math.Pi / 180
    // get my observation angles in radian
    obs := satellite.ECIToLookAngles(
        position, location, altitude,
        // get Julian date
        satellite.JDay(
            year, month, day,
            hour, minute, second,
        ),
    )
    // print satellite coordinates in angles
    fmt.Println(latlng.Latitude * math.Pi / 180)
    fmt.Println(360 + (latlng.Longitude * math.Pi / 180))
    // print my observation azimuth in angle
    fmt.Println(obs.Az * 180 / math.Pi)
}

[ThomasHabets]

Thanks. Much better. Still, I'm not getting exactly the same values as gpredict, for the same passes. Especially for elevation, which is 1-4 degrees off.

It's good enough, but strange. Looks like gpredict uses wgs72 in my version, but that doesn't seem to change the angles two decimal points down.

I assume your latitude printing is a copy-paste error? Should be * 180 / math.Pi?

Also for my sanity: Longitude, is positive West or East?

I made another test, one that doesn't have formatting errors and such that this issue had (sorry), here: https://github.com/ThomasHabets/tleservice/blob/main/test/test.go

[bclswl0827]

Thanks. Much better. Still, I'm not getting exactly the same values as gpredict, for the same passes. Especially for elevation, which is 1-4 degrees off.

It's good enough, but strange. Looks like gpredict uses wgs72 in my version, but that doesn't seem to change the angles two decimal points down.

I assume your latitude printing is a copy-paste error? Should be * 180 / math.Pi?

Also for my sanity: Longitude, is positive West or East?

I made another test, one that doesn't have formatting errors and such that this issue had (sorry), here: https://github.com/ThomasHabets/tleservice/blob/main/test/test.go

sorry for the typo, it should be * 180 / math.Pi, and I once doubted the calculation of longitude, finally, I found that the converted longitude in degrees, also need to add 360, now everything looks fine.