search

skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.39k stars 209 forks source link

Matrix Math in positionlib.ICRF.from_time_and_frame_vectors #727

Closed W-Hunter-Daboll closed 2 years ago

W-Hunter-Daboll commented 2 years ago

Lines 135 and 136 of positionlib.ICRF.from_time_and_frame_vectors call the mxv (matrix times vector) function, and (as written in the github repo) will return a matrix of nxn size.

However, a matrix multiplied by a vector should return a vector of size 1xn.

By changing line 135 from "r = mxv(RT, r)" to "r = mxv(r,RT)", this issue is resolved. Similarly, changing the order of v and RT in the velocity calculation on line 136 resolves that issue.

As an interesting side note, comparing this method to numpy.dot(r,RT) yields slightly different results - see the below example, in which all numbers are in astronomical units, being taken from the skyfield "position" class: image

I will presently be working to validate these calculated positions with a known use case.

brandon-rhodes commented 2 years ago

I will presently be working to validate these calculated positions with a known use case.

Thank you! I will wait and see what the result of your analysis is before trying anything myself. If your comparison with the known use case suggests there's an improvement to be made to Skyfield here, see if you can put together a small sample script with (a) its current output and (b) the output you would like it to have instead. Thanks!

brandon-rhodes commented 2 years ago

@W-Hunter-Daboll — Did you succeed in validating those calculated positions? I'm still interested in fixing this issue if indeed it turns out that a bad value is being computed. Thanks!

W-Hunter-Daboll commented 2 years ago

Hey @brandon-rhodes, thanks for checking back in. Unfortunately, I have not been able to undertake that validation effort, and don't expect to in the future. However, I would highly recommend investigating both the numpy.dot() discrepancy and the error in line 135 as mentioned above. Apologies that I cannot be of greater help at this time.