Geocentric altitude/az of the sun and moon.

[hidp123]

Normally, alt & az are topocentric. However, how can I get geocentric alt & az. Altitude at least, if not azimuth.

[cosinekitty]

This is an interesting question! I know an easy way to calculate the geocentric altitude, but it's not clear to me how to do the azimuth.

My approach would be to call ObserverVector to obtain a vector from the center of the Earth to the observer. You should pass false for ofdate, to select J2000 coordinates.

Then call GeoVector to obtain a vector from the center of the Earth to the celestial body in question.

Now you have two geocentric vectors, one for the observer and one for the celestial body, both in consistent J2000 equatorial coordinates.

Subtract the observer vector from the body vector to obtain a vector from the observer to the body. Now pass the geocentric observer vector and the observer-centric body vector to the AngleBetween function. This will return an angle in degrees between the two vectors, which is the body's angle with respect to the geocentric zenith. Subtract that angle from 90 degrees to obtain the geocentric altitude angle.

I hope this helps! By the way, I am curious: why do you want geocentric angles rather than topocentric angles? That is not something I expected anyone would ever want.

[cosinekitty]

OK, I will study this over the next couple of days and see if I can help you do the same thing in Astronomy Engine. Are you pretty sure that the value 7.16406779367621 (or something very close to it) is correct for your use case?

[hidp123]

It seems to be, as I posted in the screenshot above, the tables given by Yallop are only to 1 decimal place, so hard to gauge precision.
Also, JPL Horizons does not give geocentric values for az/alt, so hard to compare for verification as well.

I saw this: #53 (comment), and this: #54, not sure if it is connected to what we are doing.

[cosinekitty]

#53 and #54 are about work I did adding coordinate transforms using rotation matrices. Basically any vector in any of the supported orientation systems can be converted to a vector in any other orientation system. This could very well be of help in your case; I'm not sure yet, but certainly possible. The documentation section on coordinate transforms gives more details if you are interested.

[cosinekitty]

OK, I think I figured it out. I was calculating the vectors in J2000 coordinates (EQJ in Astronomy Engine nomenclature), but to get accurate horizontal coordinates, I need to convert them to of-date coordinates (EQD).

#!/usr/bin/env python3
import astronomy
observer = astronomy.Observer(39, -76.8)
besttime = astronomy.Time.Make(1990, 11, 19, 22, 29, 11)
geomoon = astronomy.GeoVector(astronomy.Body.Moon, besttime, True)
geosun = astronomy.GeoVector(astronomy.Body.Sun, besttime, True)
rot = astronomy.Rotation_EQJ_EQD(besttime)
rotmoon = astronomy.RotateVector(rot, geomoon)
rotsun  = astronomy.RotateVector(rot, geosun)
meq = astronomy.EquatorFromVector(rotmoon)
seq = astronomy.EquatorFromVector(rotsun)
mhor = astronomy.Horizon(besttime, observer, meq.ra, meq.dec, astronomy.Refraction.Airless)
shor = astronomy.Horizon(besttime, observer, seq.ra, seq.dec, astronomy.Refraction.Airless)
print(mhor.altitude, shor.altitude, mhor.altitude - shor.altitude)
# 7.162904852724097 -7.808661756113622 14.971566608837719

So now the discrepancy between Astronomy Engine and Swiss Ephemeris is 0.001163 degrees, or 4.2 arcseconds, which is about the amount of error I would expect.

I remembered there is an EquatorFromVector function, so we don't need the undocumented _vector2radec after all.

Also note that I changed to correcting for aberration in the GeoVector calls, because it makes a significant difference for the Sun. It makes no difference at all for the Moon, but I made them both True for consistency.

[hidp123]

Thanks for your efforts and help. I'll give this a go.

[hidp123]

I was wondering if this is similar to the work you did here #279, but instead of going higher we are going lower towards the core of the Earth. So a negative altitude of sorts, and the horizon will shrink equally as well, like a smaller sphere.

[cosinekitty]

The distinction between geocentric altitude and geodetic altitude is a different problem. This has to do with corrections for the oblateness of the Earth. If you are on a ship 45 degrees north of the equator, and you measure the altitude angle of an object above the sea level horizon, it turns out that your horizon plane is not perpendicular to a line passing through you and the geocenter. In other words, if you let a plumb line hang, and extrapolate the direction it points, that extended line will not intersect with the Earth's center. It will miss because the rotation of the Earth causes the ocean (and your sense of up and down) to bulge away from a spherical shape. Stated another way, your sense of gravitational up and down is a combination of a gravity vector and a centrifugal force vector that point in different directions.

[YingChangJ]

Hello,

I happened to see that you are still active with this problem few days ago, and try to reproduce the result with the three libraries you mentioned, they could produce pretty much the same result (astronomy-engine 0.0001 diff, and skyfield 0.00001 diff with swisseph).
As for astronomy-engine, starting with GeoVector not bad, we can get all others with rotation matrix.

import astronomy

observer = astronomy.Observer(39, -76.8)
besttime = astronomy.Time.Make(1990, 11, 19, 22, 29, 11)
geomoon = astronomy.GeoVector(astronomy.Body.Moon, besttime, True)
rot1 = astronomy.Rotation_EQJ_ECL()
rot2 = astronomy.Rotation_ECL_HOR(besttime, observer)
rot3 = astronomy.CombineRotation(rot1, rot2)
mec_hor = astronomy.RotateVector(rot3, geomoon)
print(astronomy.SphereFromVector(mec_hor))

# Spherical(lat=7.162904852724086, lon=132.47079289201685, dist=0.002715735744607922)

It is exactly the same as the one above. (while swisseph give 7.163034615148587)
(I post it here to show how the other two libraries could match, in your recent post)