skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

EarthSatellite.geocentric.subpoint Errors (possible bug) #295

Closed ykaravas closed 5 years ago

ykaravas commented 5 years ago

I have found that the geocentric.subpoint function may be off in my limited experience with it. Specifically with the subpoint.longitude value. I have written my own function to calculate this value as well which i have based off a function i have used before in c++. This function calculates the geocentric coordinates of the intersection of the vector from (0,0,0) to the satellite's geocentric location. I have used various online calculators that go from ecef to lat lon to test my geocentric coordinate output against the EarthSatellite.geocentric.position.m and they come out to be the same lat and lon so i am fairly confident my function is working properly. I have tested this for several satellites. The subpoint function, however, outputs the same latitude that online calculators output for both my intersection coordinates and EarthSatellite.geocentric.position.m coordinates but a very different longitude. It appears to be off by different values for different satellites so i dont believe it is an interpretation misunderstanding. Below is an example of a scenario i ran using the "COSMOS 2463 [+]" satellite from http://celestrak.com/NORAD/elements/musson.txt:

geocentric.subpoint calculation:

Latitude: 70deg 34' 06.9" Longitude: 71deg 01' 56.0" Elevation (m): 1002682

geocentric.position.m

array([-1790426.34037484, -1676207.77538071, 6941257.33151884])

MY projection to earth of satellite calculation in ecef:

[-1549052.735985318, -1450232.373139543, 6005482.27988889]

Lat and Lon of geocentric.position.m coordinate according to: http://www.sysense.com/products/ecef_lla_converter/index.html

70.643970, -136.887097, 1002700.12

Lat and Lon of my intersection calculation coordinate according to: http://www.sysense.com/products/ecef_lla_converter/index.html

70.660126, -136.887097, 10229.82

NOTE (altitude is different between two outputs obviously because mine is technically on the earth but even my altitude is off because i am using Earth's equatorial radius as an approximation). But as you can see, the lat lon of my ECEF coordinates that i have calculated to be the projection of the satellite onto the earth and the lat lon EarthSatellite.geocentric.position are the same while the EarthSatellite.geocentric.subpoint longitude differs.

Attached is the code I used to get the EarthSatellite.geocentric.position and subpoint as well as to calculate my projection onto earth geocentric value. (Note, the output from the first function is the input to the second function

functions.txt

ykaravas commented 5 years ago

So it turns out the problem here is the difference between ECEF and GCRS coordinate systems. I was assuming ECEF was being used as the "geocentric" term for earth satellites but that is not the case

brandon-rhodes commented 5 years ago

So it turns out the problem here is the difference between ECEF and GCRS coordinate systems. I was assuming ECEF was being used as the "geocentric" term for earth satellites but that is not the case

Thanks for the update! I'm back from a trip that prevented me from immediately responding. Do you now consider the issue resolved, in which case we can close it? Or do you still have a question or request? If so, an example script I can run would be helpful. Thanks for clarifying!

ykaravas commented 5 years ago

I would consider it resolved. However, i would make it a bit more explicit (and maybe it is and I missed it) that your geocentric values are GCRS and Not ECEF. Because in my experience, and this is not necessarily correct but unfortunately is the case, geocentric and ECEF are commonly conflated to be the same thing since ECEF is the most common geocentric coordinate system. That’s all! Thanks!

Regards, Yiannis Karavas

On Fri, Nov 1, 2019 at 23:04 Brandon Rhodes notifications@github.com wrote:

So it turns out the problem here is the difference between ECEF and GCRS coordinate systems. I was assuming ECEF was being used as the "geocentric" term for earth satellites but that is not the case

Thanks for the update! I'm back from a trip that prevented me from immediately responding. Do you now consider the issue resolved, in which case we can close it? Or do you still have a question or request? If so, an example script I can run would be helpful. Thanks for clarifying!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/295?email_source=notifications&email_token=ABR4FZAYAKIYFHHNSGQCV7DQRTU4VA5CNFSM4JEYE4WKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEC4SKFI#issuecomment-549004565, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABR4FZA3VEGADFQZQ36CLVTQRTU4VANCNFSM4JEYE4WA .

brandon-rhodes commented 5 years ago

Could you share your code so that I can see the specific calls you are making? That would help me think through how their design might improve to make the coordinate system more explicit, and would also suggest where I might improve a couple of docstrings to emphasize the same point.

ykaravas commented 5 years ago

Yea the code I used is on the issue.

On Tue, Nov 5, 2019 at 09:30 Brandon Rhodes notifications@github.com wrote:

Could you share your code so that I can see the specific calls you are making? That would help me think through how their design might improve to make the coordinate system more explicit, and would also suggest where I might improve a couple of docstrings to emphasize the same point.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/295?email_source=notifications&email_token=ABR4FZF2TTERLXGTL2F4SKTQSF7RNA5CNFSM4JEYE4WKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEDC45HI#issuecomment-549834397, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABR4FZHAM3C5FSIEHD24DLDQSF7RNANCNFSM4JEYE4WA .

brandon-rhodes commented 5 years ago

Oh, I had read past it because it had a .txt extension. I see your Python code now. It gives me

NameError: name 'sat_dict' is not defined

when I try running it, but it looks like I can address the UX issue without needing the script to run successfully.

It looks like the crux might be the line:

return geocentric.position.m

because the name position is so generic. If the reader doesn’t know that Skyfield is based fundamentally on the ICRS system of coordinates, then all they know is that a vector is being returned, but not what sort of vector.

I wonder if the name position should be deprecated in favor of a more explicit name, like gcrs_xyz or gcrf_xyz. I have always been a bit bothered by the redundancy of position.position.m anyway, and switching to an explicit name would eliminate it. (Hmm. I'll have to think through how this would work more generally, since it would mean that subclasses of ICRF would now each have a different name maybe for their coordinate x,y,z. Maybe position would remain the behind-the-scenes generic name, and things like bcrs_xyz and grcs_xyz become subclass-specific user-visible aliases?)

What do you think?

I would also have to come up, presumably, with a better name for the raw velocity.

ykaravas commented 5 years ago

Yea it wouldn't let me upload a .py file for some reason. I can't speak to the sat_dict issue you are having. I did not upload ALL the code i used so that you wouldnt have to sift through it all but rather could just see what was relevant; so i doubt it is runnable in its current state.

But yes, like is said, i think the term "geocentric" is too generic. There are many geocentric coordinate systems, the most common of which is ECEF not gcrf. So you may want to change the term "geocentric" to what you suggested. It may be worth adding an ECEF coordinate to it as well. Its very easy to convert from the latitude and longitude to an ecef position (although by that same token, a user could also just do it very easily as well). I havent done too much with satellite modeling so i dont know how useful ECEF would be but i did find myself needing it for my simulation at least.

Regards, Yiannis Karavas

On Wed, Nov 6, 2019 at 11:58 AM Brandon Rhodes notifications@github.com wrote:

Oh, I had read past it because it had a .txt extension. I see your Python code now. It gives me

NameError: name 'sat_dict' is not defined

when I try running it, but it looks like I can address the UX issue without needing the script to run successfully.

It looks like the crux might be the line:

return geocentric.position.m

because the name position is so generic. If the reader doesn’t know that Skyfield is based fundamentally on the ICRS system of coordinates, then all they know is that a vector is being returned, but not what sort of vector.

I wonder if the name position should be deprecated in favor of a more explicit name, like gcrs_xyz or gcrf_xyz. I have always been a bit bothered by the redundancy of position.position.m anyway, and switching to an explicit name would eliminate it. (Hmm. I'll have to think through how this would work more generally, since it would mean that subclasses of ICRF would now each have a different name maybe for their coordinate x,y,z. Maybe position would remain the behind-the-scenes generic name, and things like bcrs_xyz and grcs_xyz become subclass-specific user-visible aliases?)

What do you think?

I would also have to come up, presumably, with a better name for the raw velocity.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/295?email_source=notifications&email_token=ABR4FZDH47YDK2NH7W2ZXXLQSLZTXA5CNFSM4JEYE4WKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEDHHWNA#issuecomment-550402868, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABR4FZA3LZP2D54YOIAPKVDQSLZTXANCNFSM4JEYE4WA .

brandon-rhodes commented 5 years ago

There are many geocentric coordinate systems, the most common of which is ECEF…

Another difficulty is that there are many ECEF coordinate systems, including the ITRF but also including several other ECEF systems — the SGP4 module, for example, has to convert between two of them. (Just as there are several geocentric “intertial” coordinate systems, of which the GCRF is only one and geocentric J2000 is another, and Skyfield has a matrix for converting between them.)

In the code, the fact that a Geocentric coordinate is inertial and not rotational is implied by its inheriting from an inertial reference frame: Geocentric(ICRF). But, I am now dismayed to realize, the Sphinx-generated class documentation doesn’t seem to even mention superclasses in its text! It only mentions instantiation arguments:

https://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Geocentric

The summary of position types in the API docs also entirely fails to mention that they are all inertial and all aligned with the ICRS:

https://rhodesmill.org/skyfield/api.html#astronomical-positions

So I am going to take the following steps:

  1. I am going to update the summary to define them all as inertial ICRS coordinates.
  2. I am going to mention the superclass in each subclass's description.
  3. I am going to also describe each subclass as (a) intertial and (b) aligned to the ICRS axes in their descriptions.

Are there other docs that it would have helped you to have updated?

brandon-rhodes commented 5 years ago

Oh, I could also rename Geocentric to ITRF, couldn’t I? Almost no non-astronomer would know at first glance what it meant, but the docs could explain that it’s geocentric. Hmm. That would be an annoying indirection for newcomers. Maybe I would also put the word "geocentric" in the repr(), yes, I like that idea! So even someone new to references systems would see <ITRF geocentric coordinate …> and have some hope of knowing what's going on.

I'll see how the resulting code looks.

brandon-rhodes commented 5 years ago

I didn't wind up renaming any classes at this point — I tried but the result did not look good — but I made extensive updates to docs and even the repr()'s returned by classes, so hopefully with the next version of Skyfield the possible confusion will decrease.

ykaravas commented 5 years ago

Yea i figured renaming the classes might be difficult especially for those using the library as they would have to refactor too. I think additional clarification in documentation regarding the type of geocentric coordinates is adequate!

Regards, Yiannis Karavas

On Sun, Nov 10, 2019 at 10:01 PM Brandon Rhodes notifications@github.com wrote:

I didn't wind up renaming any classes at this point — I tried but the result did not look good — but I made extensive updates to docs and even the repr()'s returned by classes, so hopefully with the next version of Skyfield the possible confusion will decrease.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/295?email_source=notifications&email_token=ABR4FZC3UATHNQUOLT5PGB3QTDDIRA5CNFSM4JEYE4WKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEDVQSPA#issuecomment-552274236, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABR4FZHJ6RNBXKKW4XBBMXTQTDDIRANCNFSM4JEYE4WA .