Can Skyfield Calculate Geocentric Altitude, Azimuth and Seperation

Hi,

First of all, I been using Skyfield to complete my PhD thesis on Islamic Calendar. I able to produce a numbers of article on Icarus, Astronomy and Computing, New Astronomy and Software Impacts due to the existence of Skyfield.

Thanks Mr Brandon Rhodes for doing this.

Second, I wonder if I could calculate Altitude, Azimuth and Elongation (Seperation) from Geocentric pov. This is because there is a huge debate in South East Asia whether to use topocentric (from observer POV) or geocentric (from centre of the earth POV) in calculating lunar position for Muslim Calendar.

Hope you could enlighten me on this matter. Thanks again.

I am so glad to hear that Skyfield has been useful in your research, thank you for letting me know that it has been helpful!

I would be interested in seeing some of the debate. Is there a link where I could read some of the issues are being discussed? A real-life observer on the Earth's surface sees the Moon relative to his own position, not relative to what a hypothetical observer at the Earth's core would see, so I am wondering what the debate would involve?

hurm. This is from my point of view, as a researcher, not the official statement of the country or the organization.

The debate is rather, philosophical instead of academical. One group, called Nadlatul Ulama has long history of using moon sighting to determine Muslim calendar. Another group, called Muhammadiyah, has been using pure calculation for Muslim calendar determination, regardless of the moon visibility, sky and observer condition.

Since 1995, South East Asia countries, Brunei, Indonesia, Malaysia and Singapore has came up with a criterion for moon sighting. The criterion is a combination of moon sighting and pure calculation practice, and used to determine the Muslim Calendar. Then in 2021, they revised the criterion on based recent collected data.

The 2021 criterion has parameter of 6.4 degree of elongation and 3 degree of altitude. Whenever a new moon pass these criterion, the new Muslim month begin. Nadlatul Ulama and Muhammadiyah both agree with this parameter, however they disagree with its POV. Nadlatul Ulama, with long history of lunar crescent observation argues that the parameter is calculated with topocetric POV, since it reflect actual observation. Muhammadiyah, with tradition of pure calculation, argues that geocentric pov is better as most geocentric parameter is lower than topocentric, and easier for lunar crescent to pass the parameter.To add more confusion to this issue, there are a number of naked eye observation in 2023 to support the geocentric argument.

So, I am planning to develop a projection of Muslim calendar based on both POV, to observe whether there are a long term impact of between geocentric against topocentric computation. I could calculate it from ra dec position, just wondering is there any function from skyfield that I could use (yeah sorry, I am quite lazy haha).

Thank you for the detailed explanation! It helped me return this weekend to some of the literature on Moon observations, and I see now what is going on.

Here’s some background:

In the old days, instead of computing the relative position of the Moon and an observer with simple vector subtraction, like Skyfield does today, astronomers computed a position in two steps. First, they figured out where the Moon was from the geocenter in coordinates like ra/dec or ecliptic lat/lon; sometimes these angles could simply be read from a table in a book. Second, they applied trigonometry — in an era before calculators, every step was expensive! — to adjust those coordinate angles for the fact that a real observer is somewhere on the Earth's surface, about an Earth radius away from the geocenter in some particular direction. They called this second step "horizontal parallax" by which they meant "adjusting a coordinate for where the observer and their horizon (thus “horizontal”) really is".

Given this time-consuming two-step process, the pioneers of Moon-visibility calculations made a fateful decision. Why, they asked, should we take the second step — why should we take the effort to turn the geocentric position into a topocentric one? If we are exclusively interested in the circumstances when the young Moon is setting, when both the Sun and Moon are very near the western horizon, then the second step described above will simply push the Sun and Moon down towards the horizon a bit by a roughly constant angle for each, because the effect of the second step is simply to "lift" the observer "up" from the geocenter to the Earth's surface. So let's skip the second step, they said, and do our "curve fitting" for the binary variable "Moon visible" / "Moon not visible" with the raw geocentric coordinates, since the differences in the angles will roughly be constants.

This generates a concept that had been confusing me: "geocentric altitude". Since it's not a real physical quantity, but an artificial number created by stopping an old-fashioned calculation halfway, it took me several readings of the old literature to realize what quantity they were describing. Here, by the way, is how the skip-the-second-step technique is justified in the classic 1910 paper On the smallest visible phase of the Moon by J.K. Fotheringham:

https://adsabs.harvard.edu/full/1910MNRAS..70..527F

And here is the same quantity defined in the 1997 paper A Method for Predicting the First Sighting of the New Crescent Moon by BD Yallop:

https://astronomycenter.net/pdf/yallop_1997.pdf

So — here are my three thoughts.

First, using “geocentric altitude” isn’t the approach I would take. Since the translation from the geocenter to the Earth-surface observer changes the Sun’s position less than the Moon’s, and changes altitude more than azimuth, I would only want to attempt crescent-Moon visibility prediction with the observer’s full true topocentric position.

Second, I of course understand that you need to compute “geocentric altitude” anyway, despite its drawbacks, because the whole point of your paper will be to compare different historic and current approaches.

Third, I realized this weekend that Skyfield can quite trivially do this calculation! I kept trying to think of ways to get a topocentric location to compute the location of the geocenter instead. But there’s no need! Because Python object attributes are writeable, we can just copy over the x,y,z position from the geocenter itself and create a hybrid object that has thrown away the x,y,z that Skyfield computed for it and uses Earth’s geocenter instead:

# EDIT: THIS RESULT IS WRONG - see my further comment below

from skyfield.api import N, W, load, wgs84

ts = load.timescale()
t = ts.utc(2024, 3, 25)

planets = load('de421.bsp')
earth, moon = planets['earth'], planets['moon']
boston = earth + wgs84.latlon(42.3583 * N, 71.0636 * W)

print('Actual Moon alt/az from Boston:')

astrometric = boston.at(t).observe(moon)
alt, az, d = astrometric.apparent().altaz()

print(alt)
print(az)

print('Moon alt/az if Boston were at the geocenter:')

boston_geocenter = boston.at(t)
boston_geocenter.position = earth.at(t).position  # overwrite x,y,z coordinates
funky = boston_geocenter.observe(moon).apparent()
alt, az, d = funky.altaz()

print(alt)
print(az)

~And indeed it seems to work! The Moon is higher from the geocenter, as we would expect:~ (Edit: Yeah, these positions are way too far apart. See my further comment below.)

Actual Moon alt/az from Boston:
14deg 00' 38.7"
103deg 17' 27.2"
Moon alt/az if Boston were at the geocenter:
19deg 49' 31.7"
105deg 16' 24.9"

So, try out that code yourself, and see if it lets your project proceed. Hopefully you can determine, give some real-world sightings of the crescent Moon and also some failed sightings, whether a topocentric approach would yield even better results than the traditional shortcut approach, or whether Fotheringham was really correct way back in 1910, and the difference can be neglected.

I’ll be interested to hear about your results, if you ever think to return to this issue later and let us know how your investigation turns out!

Many thanks for looking into this issue. I will try the code you have given.

I think you have hit the nail on the head regarding how they used parallax to adjust from geocentric to topocentric. Because when I adjusted the topo altitude by the parallax it gave me very similar results to what Yallop gave in his tables in his crescent paper.

A related question now is that does Skyfield have a function where I can calculate geo and topo horizontal parallax for the moon especially, and by extension for the sun as well. And if such a function does not exist then are there any workarounds to derive this?

Thanks

A related question now is that does Skyfield have a function where I can calculate geo and topo horizontal parallax for the moon especially, and by extension for the sun as well.

Parallax is simply the difference between the two positions, so you just use the Skyfield separation_from() method.

In writing up an example, I realized that the script above gives a wrong result — the two positions are way too far apart. It turns out the trick of overwriting .position unfortunately leaves un-adjusted several secret internal values of the position object, that are used by the .apparent() method. So the result comes out wrong.

So let's try the opposite approach. Instead of generating a topocentric observer, and trying to overwrite its .position with a fake value, let's create a real and genuine geocentric position, then overwrite its public center_barycentric attribute so that it gains the ability to produce an alt/az. Thus:

from skyfield.api import N, W, load, wgs84
ts = load.timescale()
t = ts.utc(2024, 3, 25)
planets = load('de421.bsp')
earth, moon = planets['earth'], planets['moon']
boston = earth + wgs84.latlon(42.3583 * N, 71.0636 * W)

boston_observer = boston.at(t)
geocenter_observer = earth.at(t)

apparent1 = boston_observer.observe(moon).apparent()
apparent2 = geocenter_observer.observe(moon).apparent()

h_parallax = apparent2.separation_from(apparent1)
print('Horizontal parallax:', h_parallax)

# Give the `apparent2` position the Boston-oriented zenith, nadir, and
# horizon that it technically doesn't deserve.

apparent2.center_barycentric = apparent1.center_barycentric  # lie about center

print('Conventional alt / az:')
print(apparent1.altaz())

print('"Geocentric" alt / az:')
print(apparent2.altaz())

Which produces:

Horizontal parallax: 00deg 52' 21.2"
Conventional alt / az:
(<Angle 14deg 00' 38.7">, <Angle 103deg 17' 27.2">, <Distance 0.00270116 au>)
"Geocentric" alt / az:
(<Angle 14deg 52' 59.9">, <Angle 103deg 17' 38.1">, <Distance 0.00271182 au>)

Which is a much more reasonable horizontal parallax that the one given by the previous script. Try checking this new script against some other sources, and see if it gives numbers that agree.

@hidp123 Hello, it is interesting, so I tried to compare with them. With the three libraires you mentioned, the other two use earth barycenter and treats earth like a perfect sphere. To reproduce the result in swisseph, produce a frame from the geolat, and geolat manually. Or use the rotation matrix in astronomy-engine. (here rot3 is the rotation matrix from equator J2000 to altaz)

from skyfield.api import N, W, load, wgs84

ts = load.timescale()
t = ts.utc(2024, 3, 25)
planets = load('de440.bsp')
earth, moon = planets['earth'], planets['moon']
boston = earth + wgs84.latlon(1 * N, 71.0636 * W)

boston_observer = boston.at(t)
geocenter_observer = earth.at(t)

apparent1 = boston_observer.observe(moon).apparent()
apparent2 = geocenter_observer.observe(moon).apparent()

# apparent2.center_barycentric = boston_observer.center_barycentric  # lie about center
apparent2.center_barycentric._altaz_rotation = np.array(rot3.rot).T

alt, az, _ = apparent2.altaz()
print(alt.degrees, az.degrees)
# 19.85460136761465 270.1774357870765

(while swisseph is 19.85461086418853) With earth as a perfect sphere, skyfield is just 0.00001 away from swisseph. (while astronomy-engine 0.0001 away which comes from the accuracy of the moon's position)

Thus, with a geocentric lunar calendar, should temperature, pressure, obliqueness be included? (only skyfield could take obliqueness into account; swisseph and astronomy-engine could take temperature and pressure into account, skyfield could not manually add refract to non-locations; and get rid of all these make these three match)

Thanks for your interest in my work. I understand it is time consuming. Previously, in my degree year, I use Jean Meeus code geocentric and topocentric calculation. Jean Meeus uses old methodology, similarly with Yallop and Fotheringham, where he calculate geocentric parameter first and then he transform the parameter into topocentric one, using nutation, LST and so on.

I could get a geocentric value similar with Horizon using Jean Meeus algorithm, however the topocentric value is very different when I compare Jean Meeus's Topocentric with Horizon. So I am at wits ends.

My friend here Izzat Zubir (https://github.com/izzatzubir) recommend that I use negative radius for geocentric calculation. The code would be like this

For Topocentric Computation

    lokasi_cerapan = earth + Topos(latitude_degrees=lat,longitude_degrees=long, elevation_m=elevation)
    sun_astro = lokasi_cerapan.at(ts.utc((ysunset), (mosunset), (dsunset),(hsunset),(msunset),(ssunset))).observe(sun)
    sun_app_geo = sun_astro.apparent()
    sun_alt_geo, sun_az_geo, sun_distance = sun_app_geo.altaz()

For Geocentric Computation

    lokasi_cerapan = earth + Topos(latitude_degrees=lat,longitude_degrees=long, elevation_m= -6378136.6)
    sun_astro = lokasi_cerapan.at(ts.utc((ysunset), (mosunset), (dsunset),(hsunset),(msunset),(ssunset))).observe(sun)
    sun_app_topo = sun_astro.apparent()
    sun_alt_topo, sun_az_topo, sun_distance = sun_app_topo.altaz()

The value would be agreeable with most Islamic Computation Software like Accurate Time. What do you think Mr @brandon-rhodes ?

skyfielders / python-skyfield

Can Skyfield Calculate Geocentric Altitude, Azimuth and Seperation #951