New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP:adding 3rd body #381

Closed
wants to merge 18 commits into
base: master
from

Conversation

2 participants
@nikita-astronaut
Copy link
Contributor

nikita-astronaut commented May 31, 2018

Hello! I am adding the third body. To make the discussion public, I make a PR though the tests are not passing yet.

The test with the 3rd-body perturbation itself gives the wrong answer and I started to think where the mistake can be.

I wrote the check of coordinates test. In Cowell, there is the example number 12.10, where the position of Moon in xyz is given at certain time. Then this position is used to compute the perturbations.

The code that I am currently using to get the coordinates
epoch = Time(2456498.8333, format='jd', scale='tdb') moon = Orbit.from_body_ephem(Moon, epoch) moon = change_frame(moon, ICRS, GCRS)
gives the y and z components with the wrong sign. Changing these signs manually in 3rd-body tests does not make it pass, but still, the coordinates issue should be resolved...

@Juanlu001 , perhaps you have an idea?

@Juanlu001
Copy link
Member

Juanlu001 left a comment

I left some comments. General advice:

  1. Make sure the tests run in the CI as soon as possible, i.e. fix any PEP8 issues you might have. That way we can detect problems earlier in the review.
  2. Do a cleanup of unused imports (if any) and try to keep the consistency order (first stdlib, then numpy, then astropy and scipy, then poliastro)
@@ -130,3 +131,36 @@ def inertial_body_centered_to_pqw(r, v, source_body):
v_pqw = (np.array([-sin(nu), (ecc + cos(nu)), 0]) * sqrt(k / p)).T * u.km / u.s

return r_pqw, v_pqw


def change_frame(orbit, frame_orig, frame_dest):

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

I'd rather call this transform, to be consistent with Astropy method transform_to :)

----------
orbit : ~poliastro.bodies.Orbit
Orbit to transform
frame_orig : ~astropy.coordinates.BaseRADecFrame

This comment has been minimized.

This comment has been minimized.

@nikita-astronaut

nikita-astronaut Jun 10, 2018

Author Contributor

done

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

I still see it there?

orbit_dest.representation = CartesianRepresentation

return Orbit.from_vectors(orbit.attractor,
[orbit_dest.x, orbit_dest.y, orbit_dest.z] * u.km,

This comment has been minimized.

@Juanlu001

This comment has been minimized.

@nikita-astronaut

nikita-astronaut Jun 10, 2018

Author Contributor

'GCRS' object has no attribute 'xyz'

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

Sorry, it should be orbit_dest.data.xyz.


return Orbit.from_vectors(orbit.attractor,
[orbit_dest.x, orbit_dest.y, orbit_dest.z] * u.km,
[orbit_dest.v_x, orbit_dest.v_y, orbit_dest.v_z] * (u.km / u.s),

This comment has been minimized.

@nikita-astronaut

nikita-astronaut Jun 10, 2018

Author Contributor

same thing

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

In this case it's more complex, but it's still there:

In [31]: dv = GCRS(x=1 * u.km, y=0 * u.km, z=0 * u.km, v_x=0*u.km/u.s, v_y=1*u.k
    ...: m/u.s, v_z=0*u.km/u.s, representation=CartesianRepresentation, differen
    ...: tial_type=CartesianDifferential)

In [32]: dv.data.xyz
Out[32]: <Quantity [1., 0., 0.] km>

In [33]: dv.data.differentials['s'].d_xyz
Out[33]: <Quantity [0., 1., 0.] km / s>
Orbit in the new frame
"""

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

Remove extra whitespace

assert_quantity_allclose((raan * u.rad).to(u.deg) - 360 * u.deg, -0.06 * u.deg, rtol=1e-2)


def test_moon_at_right_position():

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

This test is here twice?


r, v = get_body_barycentric_posvel(body.name, epoch)
return cls.from_vectors(body.parent, r.xyz.to(u.km), v.xyz.to(u.km / u.day), epoch)
if get_velocity:

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

Instead of hardcoding the JPL ephemeris here, the user should be responsible of setting the correct ephemeris:

In [3]: Orbit.from_body_ephem(Moon)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-3-c571ef3eb3af> in <module>()
----> 1 Orbit.from_body_ephem(Moon)

~/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py in from_body_ephem(cls, body, epoch)
    164                  .format(epoch.tdb.value), TimeScaleWarning)
    165 
--> 166         r, v = get_body_barycentric_posvel(body.name, epoch)
    167         return cls.from_vectors(body.parent, r.xyz.to(u.km), v.xyz.to(u.km / u.day), epoch)
    168 

~/.miniconda36/envs/poliastro36/lib/python3.6/site-packages/astropy/coordinates/solar_system.py in get_body_barycentric_posvel(body, time, ephemeris)
    328 
    329     """
--> 330     return _get_body_barycentric_posvel(body, time, ephemeris)
    331 
    332 

~/.miniconda36/envs/poliastro36/lib/python3.6/site-packages/astropy/coordinates/solar_system.py in _get_body_barycentric_posvel(body, time, ephemeris, get_velocity)
    223             if get_velocity:
    224                 raise KeyError("the Moon's velocity cannot be calculated with "
--> 225                                "the '{0}' ephemeris.".format(ephemeris))
    226             return calc_moon(time).cartesian
    227 

KeyError: "the Moon's velocity cannot be calculated with the 'builtin' ephemeris."

In [4]: from astropy.coordinates import solar_system_ephemeris
   ...: solar_system_ephemeris.set("jpl")
   ...: 
   ...: 
Out[4]: <ScienceState solar_system_ephemeris: 'jpl'>

In [5]: Orbit.from_body_ephem(Moon)
Out[5]: 151032073 x -151032947 km x 23.5 deg orbit around Earth (♁)

so I would remove this change altogether.

Six component state vector [x, y, z, vx, vy, vz] (km, km/s).
k : float
gravitational constant, (km^3/s^2)
third_body: ~OdeSolution object

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

Why it has to be an OdeSolution object? Any callable that receives a t0 (time elapsed since start of the propagation) and returns a state vector should work, right?

third_body: ~OdeSolution object
third body that causes the perturbation
"""

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

Remove extra white line

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Jun 1, 2018

Also, this will need a rebase on current master:

screenshot-2018-6-1 poliastro poliastro

(image from https://github.com/poliastro/poliastro/network)

@Juanlu001
Copy link
Member

Juanlu001 left a comment

Detected one mistake in the Moon tests.

epoch = Time(2456498.8333, format='jd', scale='tdb')
moon = Orbit.from_body_ephem(Moon, epoch)
moon = change_frame(moon, ICRS, GCRS)
r_expected = np.array([340958.0, 137043.0, 27521.3]) * u.km

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 1, 2018

Member

These numbers are wrong, see HORIZONS:

screenshot-2018-6-1 horizons web-interface

and Curtis:

screenshot from 2018-06-01 09-52-51

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Jun 1, 2018

Added one extra comment about the Moon position tests. Notice that it should be crystal clear whether those positions are with respect to the Earth-Moon-Barycenter (EMB) or the Earth center.

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Jun 2, 2018

I made h = 0.25 in poliastro.ephem, cause otherwise testing would take too long. For 60 days, there are around 2000 time points to sample from ephem, and this takes too long. I checked, even for h = 0.5 the results are the same (even here no dependence).

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Jun 3, 2018

Hey, @Juanlu001

Now for every orbit from Curtis I check three parameters. You can see, for GEO and LEO, there is a mismatch in 1 out of 3 parameters, for HEO all the parameters fail to agree.

Do you know any other sources to validate the function?)

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Jun 4, 2018

Your code seems to be using the JPL ephemeris (via the default get_velocity=True). However, the Curtis example uses the Meeus formulas, which are available in Astropy via the builtin ephemeris. So I would not try to validate the algorithms using two different ways of computing the Moon position.

I understand that one cannot directly create an Orbit object for the Moon using the builtin ephemeris because it does not return the velocity. However, for this particular example we are only interested in the position time history, which is something that can be obtained using astropy.coordinates.get_body_barycentric and the builtin ephemeris:

http://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric.html

So I would try two things first:

  1. Rollback the get_velocity=True changes in Orbit.from_body_ephem and use get_body_barycentric directly for these tests, making sure that the builtin ephemeris is used. This way the computation of the position of the Moon will match the one done by Curtis.
  2. Run the MATLAB D.49 script to retrieve the true values instead of reading them from the plot.

Two problems might arise here:

  1. Curtis is not using DOP853, but Runge-Kutta 4-5. We might have to propagate with the same algorithm if we want true accurate validation.
  2. Astropy returns the Moon coordinates in ICRS, after stating that "Meeus algorithm gives GeocentricTrueEcliptic coordinates". However, the *TrueEcliptic frames of Astropy are no longer True in this release (see looooooong discussion in astropy/astropy#6508) and therefore there might be accuracy problems with the frame conversion there because of the nutation choice.

Lastly, as general coding advice:

  1. Rebase & force push early so you don't face conflicts later.
  2. Use pytest.mark.parametrize for the 3rd body tests, because all the code is shared and only the data is different.
@codecov

This comment has been minimized.

Copy link

codecov bot commented Jun 9, 2018

Codecov Report

Merging #381 into master will increase coverage by 0.29%.
The diff coverage is 91.3%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #381      +/-   ##
==========================================
+ Coverage   83.14%   83.43%   +0.29%     
==========================================
  Files          33       34       +1     
  Lines        1726     1757      +31     
  Branches      141      142       +1     
==========================================
+ Hits         1435     1466      +31     
  Misses        252      252              
  Partials       39       39
Impacted Files Coverage Δ
src/poliastro/ephem.py 100% <100%> (ø)
src/poliastro/twobody/perturbations.py 100% <100%> (ø) ⬆️
src/poliastro/twobody/orbit.py 95.62% <100%> (ø) ⬆️
src/poliastro/coordinates.py 87.23% <80%> (+1.86%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b659a94...a49d0a0. Read the comment docs.

@Juanlu001
Copy link
Member

Juanlu001 left a comment

Some more changes requested

----------
orbit : ~poliastro.bodies.Orbit
Orbit to transform
frame_orig : ~astropy.coordinates.BaseRADecFrame

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

I still see it there?

orbit_dest.representation = CartesianRepresentation

return Orbit.from_vectors(orbit.attractor,
[orbit_dest.x, orbit_dest.y, orbit_dest.z] * u.km,

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

Sorry, it should be orbit_dest.data.xyz.


return Orbit.from_vectors(orbit.attractor,
[orbit_dest.x, orbit_dest.y, orbit_dest.z] * u.km,
[orbit_dest.v_x, orbit_dest.v_y, orbit_dest.v_z] * (u.km / u.s),

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

In this case it's more complex, but it's still there:

In [31]: dv = GCRS(x=1 * u.km, y=0 * u.km, z=0 * u.km, v_x=0*u.km/u.s, v_y=1*u.k
    ...: m/u.s, v_z=0*u.km/u.s, representation=CartesianRepresentation, differen
    ...: tial_type=CartesianDifferential)

In [32]: dv.data.xyz
Out[32]: <Quantity [1., 0., 0.] km>

In [33]: dv.data.differentials['s'].d_xyz
Out[33]: <Quantity [0., 1., 0.] km / s>
@@ -40,11 +40,13 @@ def test_compute_soi():
assert_quantity_allclose(r_SOI, expected_r_SOI, rtol=1e-1)


'''
@pytest.mark.parametrize("missing_body", [Moon, Pluto])

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

Why did you have to comment out this test? Also, in this case you can use pytest.mark.xfail or pytest.mark.skip https://docs.pytest.org/en/latest/skipping.html#skipping-test-functions

This comment has been minimized.

@nikita-astronaut

nikita-astronaut Jun 10, 2018

Author Contributor

I couldn't figure out the meaning of this test and wanted to discuss it with you. after my changes it does not pass anymore

from astropy.tests.helper import assert_quantity_allclose
from poliastro.twobody import Orbit
# from poliastro.coordinates import transform

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

Can we remove this commented out import?

for ti, ri, vi in zip(np.linspace(0, tof, 400), r, v):
_, _, inc, raan, argp, _ = rv2coe(Earth.k.to(u.km**3 / u.s**2).value, ri, vi)

# fighting %2pi

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

I guess this is done to put argp, raan and inc in the range [-pi, pi]? In that case, you could use http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html

In [45]: a = Angle([-20.0, 150.0, 350.0] * u.deg)

In [46]: a.wrap_at(180 * u.deg)
Out[46]: <Angle [-20., 150., -10.] deg>

In [47]: import numpy as np

In [48]: a.wrap_at(np.pi * u.rad)
Out[48]: <Angle [-20., 150., -10.] deg>
incs.append((inc + np.pi) % (2 * np.pi) - np.pi)

# averaging in order to get agreement with Curtis
inc_f, raan_f, argp_f = np.mean(incs[-5:]), np.mean(raans[-5:]), np.mean(argps[-5:])

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

Why do you use the last 5 values only? This deserves a comment

@@ -162,9 +162,12 @@ def from_body_ephem(cls, body, epoch=None):
warn("Input time was converted to scale='tdb' with value "
"{}. Use Time(..., scale='tdb') instead."
.format(epoch.tdb.value), TimeScaleWarning)
if get_velocity:

This comment has been minimized.

@Juanlu001

Juanlu001 Jun 10, 2018

Member

I still would remove this functionality, and make the user responsible of setting the correct ephemerides from the outside.

nikita-astronaut added some commits Jun 10, 2018

@Juanlu001
Copy link
Member

Juanlu001 left a comment

Looks good to me! I will merge this manually to squash some commits and cleanup history.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Jun 11, 2018

Merged manually :) Thanks a lot @nikita-astronaut! This was a tough one 💪

@Juanlu001 Juanlu001 closed this Jun 11, 2018

@nikita-astronaut nikita-astronaut deleted the nikita-astronaut:3rd_body branch Jun 12, 2018

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Aug 7, 2018

Before merge, we had to make a very important comment!

The reason why tests disagreed was than in Curtis, they claimed that argp_initial is given in deg, but actually in their simulations they forgot to transform and took value on rad, which caused huge discrepancy. As we looked at their matlab code and used the value of argp in rad, everything agreed perfectly.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment