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

added draft for mean_motion propagator #322

Merged
merged 11 commits into from Feb 27, 2018

Conversation

2 participants
@nikita-astronaut
Copy link
Contributor

nikita-astronaut commented Feb 22, 2018

Hello @Juanlu001

There is a lot to correct, I believe, but there's a mean_motion function. Now iss and halleys[0] propagate without falling and results of propagation agree with the kepler propagator where is doesn't fall.

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 22, 2018

I must admit I haven't tested the case ecc == 1. How can I do this? How to import such an orbit?

1 similar comment
@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 22, 2018

I must admit I haven't tested the case ecc == 1. How can I do this? How to import such an orbit?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 22, 2018

Now iss and halleys[0] propagate without falling and results of propagation agree with the kepler propagator where is doesn't fail.

Wow, impressive! I will try to take a closer look this weekend.

For the ecc == 1 case, you have a parabolic classmethod:

http://docs.poliastro.space/en/latest/api.html#poliastro.twobody.orbit.Orbit.parabolic

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 22, 2018

The code falls, not all the tests can be passed. I faced the problem that nu_to_M function returns M in radians or dimensionless. I don't understand, why so. Will try to find out.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

The code falls, not all the tests can be passed.

The code fails because of PEP8 errors, we still don't know whether the tests are fine. Perhaps we should split that part as well.

I faced the problem that nu_to_M function returns M in radians or dimensionless.

You will probably need to do something like this:

with u.set_enabled_equivalencies(u.dimensionless_angles()):

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

This needs a rebase on current master.

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 23, 2018

Hey, @Juanlu001

Now the code fails because it doesn't pass the test_sample_angle_zero_returns_same() test. The reason is the following: now the propagate functions transfers nu to M and back. Because of these conversions, the error occurs. For instance, the x coordinate is in one case 1131.3399999999997 (after propagation), in the other 1131.34 (before propagation). And the assertion fails. How do you think it would be best to fix it?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

I see! I think the best thing is to use assert_quantity_allclose:

In [1]: from astropy.tests.helper import assert_quantity_allclose

In [2]: from astropy import units as u

In [3]: assert_quantity_allclose(1 * u.km, 1000 * u.m)

In [4]: assert_quantity_allclose([1] * u.km, [1000] * u.m)

In [5]: assert_quantity_allclose([1, 2] * u.km, [1000, 2000] * u.m)
@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 23, 2018

Ok! So you suggest me to also change tests files?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

Of course I will examine those changes closely, but in this case I think it's necessary.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

Also, I will require you to add tests that fail with the current propagator and that succeed with the propagator that you're adding. You can use the Halley and ISS orbits as you already saw.

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 23, 2018

@Juanlu001 , hey! I faced the problem I cant really explain.
In function test_sample_angle_zero_returns_same() I make (for debug)
print(np.allclose([(ss0.r[0] / u.km).value], [(rr[0].get_xyz()[0] / u.km).value])) assert np.testing.assert_allclose([(ss0.r[0] / u.km).value], [(rr[0].get_xyz()[0] / u.km).value])

Print returns True, but assert fails. Can it be some trick of dimension wrappers and checkers that are everywhere?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 23, 2018

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Feb 23, 2018

It uses the numpy function assert_allclose. But the behavior is the same, both asserts fail.
print(a) print(b) assert_quantity_allclose(a, b)
gives
6672.423 km 6672.423 km Traceback (most recent call last): File "./src/poliastro/tests/tests_twobody/test_sample.py", line 98, in <module> test_sample_angle_zero_returns_same() File "./src/poliastro/tests/tests_twobody/test_sample.py", line 26, in test_sample_angle_zero_returns_same assert assert_quantity_allclose(a, b) AssertionError

@codecov

This comment has been minimized.

Copy link

codecov bot commented Feb 24, 2018

Codecov Report

Merging #322 into master will increase coverage by 0.35%.
The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #322      +/-   ##
==========================================
+ Coverage    80.2%   80.56%   +0.35%     
==========================================
  Files          30       30              
  Lines        1440     1461      +21     
  Branches      113      114       +1     
==========================================
+ Hits         1155     1177      +22     
  Misses        252      252              
+ Partials       33       32       -1
Impacted Files Coverage Δ
src/poliastro/twobody/propagation.py 90.36% <100%> (-0.12%) ⬇️
src/poliastro/twobody/orbit.py 95.34% <100%> (+2.38%) ⬆️

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 0492f60...c4a5cf1. Read the comment docs.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 24, 2018

Fixed :) When you add some tests and all checks pass, I will do a complete review of the code so you can make some modifications as needed, and then this is ready to merge.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 24, 2018

@Juanlu001
Copy link
Member

Juanlu001 left a comment

I added some comments but I like this very much. Looking forward to have this in poliastro, it's a very important change.

@@ -24,7 +26,7 @@ def test_propagation():

ss0 = Orbit.from_vectors(Earth, r0, v0)
tof = 40 * u.min
ss1 = ss0.propagate(tof)
ss1 = ss0.propagate(tof, method=method)

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Awesome!

@@ -6,7 +6,8 @@

from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import propagate
from poliastro.twobody.propagation import propagate, kepler, mean_motion

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Is propagate still needed here?

@@ -6,7 +6,8 @@

from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import propagate
from poliastro.twobody.propagation import propagate, kepler, mean_motion
import numpy as np

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Same for this

@pytest.mark.parametrize("time_of_flight,function", [
(1 * u.min, propagate),
(40 * u.min, propagate),
@pytest.mark.parametrize("time_of_flight,method", [

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Prefer stacking @pytest.mark.parametrize here to get all the combinations, as explained here: https://docs.pytest.org/en/latest/parametrize.html#pytest-mark-parametrize-parametrizing-test-functions (bottom of the section)

@pytest.mark.parametrize("time_of_flight,function", [
(6 * u.h, propagate),
(2 * u.day, propagate),
@pytest.mark.parametrize("time_of_flight,method", [

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Same as above

nu = M_to_nu(M, ecc)

# parabolic orbit
else:

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

I think this is a good addition, but needs testing.

# get the initial mean anomaly
M0 = nu_to_M(nu0, ecc)
# elliptic or hyperbolic orbits
if ecc != 1.0:

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Just for curiosity: what is the asymptotic behavior of this when the eccentricity is very close to 1.0? We should add a # TODO: here to remind us that this needs better validation. Perhaps in the original paper there's some information.

@@ -10,12 +11,13 @@
from poliastro.constants import J2000
from poliastro.bodies import Sun, Earth
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import cowell
from poliastro.twobody.propagation import cowell, kepler, mean_motion

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Also, can you add a test that demonstrates that the new mean_motion propagator works with the Halley orbit, and another one marked as @pytest.mark.xfail with the kepler propagator?

# get the initial mean anomaly
M0 = nu_to_M(nu0, ecc)
# elliptic or hyperbolic orbits
if not np.isclose(ecc, 1.0, rtol=1e-06):

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Can you explain why the parabolic test failed before this line was changed?

@@ -65,6 +68,23 @@ def test_propagation_hyperbolic():
assert_quantity_allclose(norm(v), expected_v_norm, rtol=1e-3)


def test_propagation_mean_motion_parabolic():
# example from Howard Curtis, section 3.5, problem 3.15

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Please specify 3rd edition

@@ -65,6 +68,23 @@ def test_propagation_hyperbolic():
assert_quantity_allclose(norm(v), expected_v_norm, rtol=1e-3)


def test_propagation_mean_motion_parabolic():
# example from Howard Curtis, section 3.5, problem 3.15
p = 13200.0 * u.km

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Can you explicitly write 2 * 6600 here?

p = 13200.0 * u.km
_a = 0.0 * u.deg
orbit = Orbit.parabolic(Earth, p, _a, _a, _a, _a)
orbit = orbit.propagate(0.8897 * 3600.0 / 2.0 * u.s, method=mean_motion)

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

What about 0.8897 * u.h / 2? :)

assert_quantity_allclose(nu0, np.deg2rad(90.0), rtol=1e-4)

orbit = Orbit.parabolic(Earth, p, _a, _a, _a, _a)
orbit = orbit.propagate(36.0 * 3600.0 * u.s, method=mean_motion)

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

Same as above, use 36 * u.h


orbit = Orbit.parabolic(Earth, p, _a, _a, _a, _a)
orbit = orbit.propagate(36.0 * 3600.0 * u.s, method=mean_motion)
assert_quantity_allclose(np.sqrt(np.sum(np.array(orbit.r.to(u.km).value) ** 2)), 304700.0, rtol=1e-4)

This comment has been minimized.

@Juanlu001

Juanlu001 Feb 27, 2018

Member

You can use from poliastro.util import norm here

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Feb 27, 2018

This is good to go. Thanks a lot for this pull request!

@Juanlu001 Juanlu001 merged commit 60d463f into poliastro:master Feb 27, 2018

5 checks passed

codeclimate All good!
Details
codecov/patch 100% of diff hit (target 80.2%)
Details
codecov/project 80.56% (+0.35%) compared to 0492f60
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@wafflebot wafflebot bot removed the 2 - In Progress label Feb 27, 2018

@nikita-astronaut nikita-astronaut deleted the nikita-astronaut:iss265 branch Feb 27, 2018

@Juanlu001 Juanlu001 referenced this pull request Aug 6, 2018

Merged

writing my GSoC'18 report #17

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