Elevation above the horizon?

[Krytic]

Hey all!

I'm using Skyfield to predict the positions of some satellites I'm interested in. I use the following code snippet to detect if a satellite will be over the horizon at time t:

    me = Topos(...)
    difference = satellite - me

    topocentric = difference.at(t)

    alt, az, distance = topocentric.altaz()

    if alt.degrees >= 0:
        pos = satellite.at(trange) # trange is a vector of times greater than t 
        path = pos.subpoint() # type(path) = Topos

What I would like to do now is to calculate the elevation as an angle above the horizon as seen by me at each point. I've been using geometry to derive an expression for elevation, but this has been fruitless when comparing to known elevations from software like gpredict.

Does Skyfield have some sort of inbuilt method to do this?

[brandon-rhodes]

Could you supply a little script that I can run locally to compute the position of a satellite you are interested in, and then name exactly the angle you want? "For this satellite at this time, I want it to print 63°" would be specific enough that I could give some advice. I can't quite tell what "elevation as an angle above the horizon" means without a specific example, because it sounds like the definition of the "altitude" number you already have?

[Krytic]

Here's an example for the satellite ASTROCAST 0.1 at 11:06PM (UTC) April 7, 2019, for the station at my University at 36.8534 S, 174.7684 E.

Gpredict has this plot for it's az/el . 11:06PM UTC is the elevation in the center, which is ~22 degrees. In the sample code below, get_elevation_angle(pos) should return about 22.

This is a MWE of the script I'm currently using.

from skyfield.api import Topos, load, utc
import datetime
import numpy as np

min_range = 15
now = datetime.datetime.utcnow()
mtime = now + datetime.timedelta(minutes=min_range)

sat = 'ASTROCAST 0.1'

### This is the function I am aiming to implement.
def get_elevation_angle(pos):
    pass

ts = load.timescale()
t = ts.utc(mtime.replace(tzinfo=utc))
mins = np.arange(1, min_range+5)

trange = ts.utc(now.year, now.month, now.day, now.hour, mins, now.second)

stations_url = 'http://celestrak.com/NORAD/elements/active.txt'
satellites = load.tle(stations_url)
satellite = satellites[sat]

auckland = Topos('36.852664 S', '174.768017 E')
difference = satellite - auckland

topocentric = difference.at(t)

alt, az, distance = topocentric.altaz()

if alt.degrees >= 0:
    # satellite is above the horizon
    pos = satellite.at(trange)

    get_elevation_angle(pos)

The aim is essentially to recreate the plot GPredict uses, to exploit for further functionality, albeit discretized.

[brandon-rhodes]

If I set t to the time you mentioned:

t = ts.utc(2019, 4, 7, 23, 6)

and then ask for the altitude:

print(alt.degrees)

I get 21.424982961350132 which is right around the 22° you see on the chart. Could you double-check and see if you get a similar result?

[Krytic]

Yes, I do, but that wasn't quite what I was after. I was computing positions using a range of times. In the MWE I gave, I wanted the elevation of pos. Since pos is an instance of Topos, I couldn't do anything like using altaz() or pos.altitude since I was after an angle.

What I ended up doing is, in place of get_elevation_angle(pos):

for tp in trange:
    diff = satellite - auckland
    diff = diff.at(tp)

It's not the cleanest solution, but it works.