Welcome to Skyfield discussions! (To ask a question, click "New Discussion" rather than add a comment to this topic)

[brandon-rhodes]

Welcome!

I have activated this GitHub “Discussions” tab to make Skyfield more welcoming for contributors and users!

• Ask questions about exceptions in your code, if you are not sure whether you are looking at a bug in Skyfield or a problem with some other cause.
• Ask questions that involve a mix of astronomy, Python, and general personal confusion — the kinds of questions that, alas, tend to be closed as too broad or open-ended (or as lacking enough research) on popular Q&A sites for programming.
• If you need a new astronomy feature in Python but are not sure whether it’s a good fit for Skyfield or not, feel free to describe the problem you are trying to solve and see what feedback the community and maintainers might offer.

[davidmikolas]

Skyfield rocks, I'm all in!

[davidmikolas]

There is an interesting Astronomy Stack Exchange question that references Skyfield and though I think https://rhodesmill.org/skyfield/searches.html can help I'm not sure how. I'm also not sure if posting interesting SE questions here is appropriate or not, let's see what happens...

https://astronomy.stackexchange.com/questions/41482/how-to-find-the-local-azimuth-of-the-highest-point-of-the-ecliptic

[brandon-rhodes]

Well, that’s a bit dispiriting — I already gave an initial answer to that, but on Stack Overflow; it looks like the question is being posted multiple places without any thought to the duplication of effort that might entail.

https://stackoverflow.com/questions/66232657/how-to-find-the-local-horizon-longitude-of-the-highest-point-of-the-ecliptic

I've pasted my answer on the astronomy site as well, though I’m not sure if doing the work of pasting an answer everywhere a duplicate question is asked is really the appropriate response ­— and I am a bit surprised that the Stack Exchange system doesn’t notice duplicate questions across its sites.

[davidmikolas]

@brandon-rhodes I'll go after the cross-posts for you, those are always discouraged and n-1 of them must be closed to prevent answer fragmentation. Thank you for you speedy response!

[smp55]

Love Skyfield. Thanks so much! I know it just barely came out, but I was wondering if you have started looking into SGP4-XP? (I know you also maintain the SGP4 library, and that GitHub doesn't have a discussion page)

[brandon-rhodes]

I'm glad you're enjoying Skyfield!

I guess I should go turn on GitHub discussions for the SGP4 library, since it's working fine here. There, it’s turned on!

I hadn't heard about SGP4-XP. I can't see the zipfile referenced from discussions online, which apparently used to be here:

https://www.space-track.org/documentation#/sgp4

If rumor is correct that only a compiled version is supplied, then I imagine someone will want to write a little Python library that wraps it up for use. I suspect it will be cleanest if that library is separate from the old sgp4 since they by definition won't use each other's code.

[smp55]

That's right. Some folks at my company developed it on behalf of the government (now Space Force, I believe), which is how I came to know about it. It's still located at that link (you have to create a space-track user account to access that page, otherwise it redirects to the API page). It still says SGP4, but the new -XP stuff is in there as well (v8.0). Only the DLLs are available currently (not sure what the plan for the source is). Unfortunately, that's a bit beyond my current Python skill level, so hopefully somebody else is on it!
FYI, the plan is that space-track will start publishing TLEs in both formats (SGP4 and SGP4-XP), differentiating using column 63 in line 1: Zero or blank for SGP4, now 4 for SGP4-XP. As I understand it, it's supposed to increase accuracy by at least an order of magnitude everywhere from LEO out to cis-lunar.

[brandon-rhodes]

Oh, wow, so you're at the same place it was developed! Neat. I know that the inaccuracy of SGP4 frustrates users, so it's great news that a better library is on the way. Please comment again if you happen to hear about source being available someday in the future.

[smp55]

One more question, since I'm here: is there a way to calculate az/el rates in skyfield? I didn't see it, though I'm certainly no expert in the field. Right now I have to use astropy for that. Nothing against astropy, but I can do everything I need with skyfield, that would be better!

[brandon-rhodes]

That's a good question! There are two obstacles right now to getting az/el rates:

1. There's no method to do it. We should add one. Could you add an issue? There's a method frame_latlon() that gives angles; and a method frame_xyz_and_velocity() that does velocity; but nothing like “frame latlon and velocity” that would do both.
2. I don't know the math. I'll need to look up some spherical trigonometry or else use a pencil and paper to remember the right approach.

Here's what I just tried:

from skyfield import api
from skyfield.constants import tau
from skyfield.functions import to_spherical

ts = api.load.timescale()
t = ts.utc(2021, 1, 28, 22, 22)

e = api.load('de421.bsp')
geographic = api.wgs84.latlon(34.0232, -84.3616)

a = (e['earth'] + geographic).at(t).observe(e['mars'])
print(a)
print(a.position.km)
print(a.velocity.km_per_s)

r, v = a.frame_xyz_and_velocity(geographic)

distance, alt, az = to_spherical(r.km)
radial_velocity, alt_velocity, az_velocity = to_spherical(v.km_per_s)

print(distance, 'km')
print(alt / tau * 360, 'degrees altitude')
print(az / tau * 360, 'degrees azimuth')
print(radial_velocity, 'radial velocity km/s')
print(alt_velocity / tau * 360, 'altitude degrees/s [WRONG]')
print(az_velocity / tau * 360, 'azimuth degrees/s [WRONG]')

And the result:

<Astrometric ICRS position and velocity at date t center=WGS84 latitude +34.0232 N longitude -84.3616 E elevation 0.0 m target=499>
[1.31926881e+08 1.02337698e+08 4.85233963e+07]
[ 0.81445906 22.30196025 10.58144965]
173874167.6684296 km
65.66748916731063 degrees altitude
132.1614247541092 degrees azimuth
24.698337015907658 radial velocity km/s
30.786040338956006 altitude degrees/s [WRONG]
78.48699757379347 azimuth degrees/s [WRONG]

Obviously turning the velocity vector into angles is the wrong approach, since it gives nonsense numbers; the angular velocities will need to be affected somehow by the position itself, since the same velocity in the sky has different names depending on where its origin falls on the celestial sphere.

Another problem is that the above code works only for astrometric positions. I don't even know the formulae that would apply the apparent-position effects, refraction and aberration, to a velocity vector, but maybe they exist. I'll check the Explanatory Supplement to the Astronomical Almanac in a week or two when I next visit my office.

[smp55]

If it helps, the process my co-worker found using astropy does the following:

1. Create Skyfield EarthSatellite object
2. Generate position and velocity vectors using ITRF_position_velocity_error method
3. Convert velocity vector using astropy.coordinates.CartesianDifferential
4. Generate cartesian pos/vel vectors using astropy.coordinates.CartesianRepresentation (including the differentials vector)
5. Create ITRS object using astropy.coordinates.ITRS and the cartesian vectors
6. Transform to alt/az frame using transform_to method and astropy.coordinates.AltAz with the observer's location
7. Alt/Az can be had with alt and az.degree methods
8. Alt and Az rates are found using the pm_alt and pm_az_cosalt methods, converting to degree per second units

[brandon-rhodes]

@smp55 — I was able to find time this week to work on a routine to compute the rates of change of azimuth and altitude! Here's example code for using the new feature:

from skyfield import api

ts = api.load.timescale()
t = ts.utc(2021, 1, 28, 22, 22)
e = api.load('de421.bsp')
geographic = api.wgs84.latlon(34.0232, -84.3616)

a = (e['earth'] + geographic).at(t).observe(e['mars'])

alt, az, distance, alt_rate, az_rate, range_rate = (
    a.frame_latlon_and_rates(geographic))

print(distance, 'km')
print(alt.degrees, 'degrees altitude')
print(az.degrees, 'degrees azimuth')
print(alt_rate.degrees.per_minute, 'degrees per minute')
print(az_rate.degrees.per_minute, 'degrees per minute')
print(range_rate.km_per_s, 'kilometers per second')

You will need the development version of Skyfield installed to try the new feature out:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

If you get the chance to try it out, let me know what you think of the results! I'll note that HORIZONS usually multiplies the azimuth rate by the cosine of the altitude, so that it's a measure of arcseconds across the surface of the sky, rather than a measure of how fast the actual azimuth number is changing (which changes much faster towards the zenith).

[smp55]

Thanks so much for this! I finally got around to updating and using it in my code, and it works perfectly. I was a little confused, because the method name is lat/lon, when really it's alt/az, but the math works perfectly!

[flight404]

Hello! I just stumbled across Skyfield while looking for an alternative to pymeeus. It seems great, but I still have much to learn. My questions will probably feel a little naive so be gentle :)

Background: I am working on a world globe data-vis setup in Houdini (showing all manner of GPS-based data). I wanted to add the feature of being able to have the correct sun and moon positions based on the time stamp which is how I found my way to Skyfield.

The good news is that I was able to get it working and I put in the date of the 2024 April 8th solar eclipse and I have a beautiful accurate result.

However, plugging in dates of future eclipses in 2028 show a correct sun-moon-earth result but the country alignment is off. This eclipse should be hitting central Australia.

This leads me to the conclusion that there is a lot about coordinate frames and ecliptic vs equatorial vs galactic coordinates that elude me.

Question: I am currently rotating the earth around its axis (in my case, the y-axis) 1 full rotation every 24 hours. This seems like what I should be doing, but maybe my reasoning is incorrect? I tried seeing if there was a specific offset or rotation multiplier that would fix my orientation issue for a few known eclipses but I didnt notice any specific pattern in how off my Earth rotation was.

Any tips or thoughts would be appreciated. Thanks!

[brandon-rhodes]

@flight404 — I think your approach is hitting the problem that 24 hours brings the Earth back to facing the Sun again, but that every day of the year the Sun is another 1/365th of the way around us in space (from our point of view), so that every 24 hours we're facing a slightly different direction. Take a look at t.gmst and see if the result, which is in hours 0–24, gives you an Earth orientation that ignores the Sun and just shows you which way the Earth is facing in space!