skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.39k stars 209 forks source link

Finding max elevation of LEO Satellite from orbital inclination and altitude only [No TLEs Allowed] #789

Closed amanchokshi closed 1 year ago

amanchokshi commented 1 year ago

Hey Brandon,

I know this is a long shot, but have been a fan of skyfield for ages, and thought that it may be worth asking here. I'm attempting to filter a whole bunch of satellites from a space-track.org satellite catalog, which has orbital inclination, apogee and perigee. I'm hoping to determine which of these satellites are visible at the Geographic South Pole, and what their max elevation would be, without downloading TLEs for all thousands of them.

More on the project: I'm spending a year at the South Pole, working with the South Pole Telescope, and in my free time I'm attempting to determine all satellites which pass through our observing field, and determine whether they're visible in our data.

I know how to use skyfield with TLEs and brute force this problem, but I think there must be a more elegant solution. Any ideas would be appreciated!

Cheers, Aman

brandon-rhodes commented 1 year ago

I'm glad you've enjoyed using Skyfield—and, wow, a year at the South Pole will be quite an adventure! My favorite book as a child was Surface at the Pole about the exploration of the Arctic, but that involved missions of only a few weeks each, and thus only glimpses of the polar summer and winter. You'll get to see each of them at full length!

I know of no easy way to categorize satellites without their TLEs, so my guess would be that you should go ahead and plan on downloading the full TLE catalog?

You might want to look at the SGP4 library that Skyfield uses under the hood:

https://github.com/brandon-rhodes/python-sgp4

If you just load up TLEs one after another, it should be quite fast. Once a TLE is loaded, you should find some ‘mean elements’ giving an idea of the shape of its orbit, if its orbit were a Kepler ellipse—look for the subsection ‘Mean Elements From Most Recent Propagation’ in the Attributes section of the docs:

https://pypi.org/project/sgp4/#attributes

I wonder if the rough shape of the orbit plus some basic geometry could narrow down which satellites needed further study as ones that could pass within view of the telescope?

brandon-rhodes commented 1 year ago

Re-reading the title of your question: are you asking, not how to get the two parameters, but for the geometry formula for maximum height from the pole given those two parameters? I'd be tempted to start by considering the planet to be the unit circle in the conventional xy plane, and the antarctic base to have the coordinates (0,-1), putting it down at the circle's bottom. Given orbital inclination i and altitude a and arbitrarily imagining the orbit slashing from upper-left to lower-right through the origin, then I think the xy coordinates of the satellite at maximum position south would be:

$x = (1 + a) \sin i$ $y = (1 + a) \cos i$

You could then make those coordinates relative to the South Pole by adding 1 to y, I think? Then take the arctangent to compute an angle from the horizon?

amanchokshi commented 1 year ago

Yeah! I was asking how to determine the maximum elevation (angle above horizon), given an orbital inclination and and altitude of orbit above the earth.

I've worked out the trig for a circle with radius 6358 (Earth radius at Pole), and an inclined orbit. The function seems to return sensible results. I only asked here in case you could think of an elegant skyfield way of determining the max elevation.

def calc_sat_params(alt,i):
    """
    Calcualte satellite Elevation, Distance.

    Parameters
    ----------
    alt : float
        Satellite Altitude above the Earth's surface [km]
    i : inclination
        Inclination of Satellite orbit [deg]

    Returns
    -------
    dist, elevation
    """
    r_earth = 6358. ## km, at pole

    alt_horiz =r_earth*(np.cos(np.deg2rad(abs(90-i)))-1)
    horiz_dist = r_earth*np.tan(np.deg2rad(abs(90-i)))

    dist_to_sat = np.sqrt(horiz_dist**2+(alt-alt_horiz)**2)
    satel = np.rad2deg(np.arcsin((alt-alt_horiz)/dist_to_sat))

    return dist_to_sat,satel

Thanks for the suggestions on mean elements, I hadn't even considered it!

brandon-rhodes commented 1 year ago

Ah, thanks for clarifying the question—no, Skyfield doesn't have any helpers I can think of that would be better here than good old-fashioned trig functions applied to the geometry of your situation.

Note that using the max altitude might include satellites that don't in fact crest the horizon high enough, if they have deeply elliptical orbits and the max altitude happens way up on the north end of the orbit rather than the south end.

Are the planet's ring of geostationary satellites visible from the South Pole thanks to refraction at the horizon? I don't have an intuition of whether refraction would make up for the fact that they lie an earth-radius north of the South Pole.

I'm going to close this issue since the actual answer is "nope, Skyfield won't be any help", but welcome any further comments, and would be very interested to hear any future updates about how the project goes. Wow, you're going to the South Pole!