Handle polar positions in height calculation

[jamesgrimmett]

I found an issue with the calculation of heights for polar positions while working on #839. The source of the problem is that the formula used to calculate height (height_au = R / cos(lat) - aC) is not valid when lat=+/-90, resulting in incorrect elevation in the output, e.g.,

from skyfield.api import load, wgs84, EarthSatellite

ts = load.timescale()
t = ts.utc(2023, 3, 2, 4, 13, 0)

pos = wgs84.latlon(90.0, 0.0, 0.0).at(t)

wgs84.geographic_position_of(pos)
<GeographicPosition WGS84 latitude +90.0000 N longitude 69.2181 E elevation -165111.6 m>

I've added a check for polar positions during height calculations, and in those calculations use height_au = abs(xyz_au[2]) - self.semi_minor_axis.au instead, which I think resolves the issue;

wgs84.geographic_position_of(pos)
<GeographicPosition WGS84 latitude +90.0000 N longitude 69.2181 E elevation -1.0e-09 m>

[brandon-rhodes]

That was a fun problem! It's always a useful exercise when a bug gets me to pull out a pencil and paper and work out the details of formulae that I had pasted into Skyfield in haste.

In this case the formulae came from:

https://celestrak.org/columns/v02n03/

Dr. Kelso's expression for h, down at the very bottom of the article, doubtless worked flawlessly in his context, as I don't imagine an Earth satellite ever crosses so exactly over the pole as to expose the formula’s weakness with latitudes within a few hairs of a fraction of 90°.

Thank you for contributing a possible solution! But I wanted to see whether there was a way to solve the imprecision without creating a special case or needing two paths through the code:

https://github.com/skyfielders/python-skyfield/commit/a355465779a2a45eeda8ae5b4d942271be9de3e5

Happily, I found that inside of the arctan2() call of the _compute_latitude routine were the length of both legs of the big triangle whose three corners are the target's location (x,y,z), the point on the Earth's axis (0,0,z), and the center of the radius of curvature way down the Earth's axis on the opposite side of the equatorial plane. So instead of trying to compute the length of the hypotenuse by using cos(lat) and one leg, which becomes very imprecise near 90°, we can use both legs with Pythagoras instead, and get a very high precision result in all cases.

It requires returning an additional value from _compute_latitude(), which incurs a slight cost, but it's worth the far better results that the additional value makes possible.

Thanks so much for reporting this, and for showing one way the problem could be solved—it was very helpful to start from a working solution!

I'll plan to release a new Skyfield version by the end of the month; maybe sooner. In the meantime, here's how to install the development version to test out the fix on your own box if you'd like:

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

[jamesgrimmett]

Excellent!

Thanks for sharing your working and thought process, it's very useful to see. I also felt that my proposed solution was a little less than ideal, so I'm glad that you found a tidier solution.

I have been referring to Dr. Kelso's columns quite a bit while learning, so I will keep this alternative solution in mind (it's also a good reminder that there is always a chance the existing formulae can be improved or reworked!).

I'll rebase the branch for #839 to include this new work. Thanks!