@@ -0,0 +1,1415 @@
+import numpy as np%0A%0Afrom astropy import units as u%0A%0Afrom poliastro.bodies import Earth%0Afrom poliastro.twobody import State%0Afrom poliastro.twobody.propagation import cowell%0Afrom poliastro.util import norm%0A%0A%0Adef test_inclination():%0A    a_0 = 7000.0%0A    i_0 = (28.5 * u.deg).to(u.rad).value%0A    f = 3.5e-7%0A%0A    k = Earth.k.decompose(%5Bu.km, u.s%5D).value%0A%0A    def a_d(t0, u_, _):%0A        r, v = u_%5B:3%5D, u_%5B3:%5D%0A%0A        beta = np.pi / 2 * np.sign(r%5B0%5D)  # Change with out-of-plane velocity%0A        #beta = np.pi / 2 * np.sign(r%5B1%5D)  # Change at node crossing%0A%0A        # DEBUG%0A        ss = State.from_vectors(Earth, r * u.km, v * u.km / u.s)%0A        print(beta, ss.inc.to(%22deg%22))%0A        # END DEBUG%0A%0A        w_ = np.cross(r, v) / norm(np.cross(r, v))%0A        accel_v = f * np.sin(beta) * w_%0A        return accel_v%0A%0A    # Retrieve r and v from initial orbit%0A    s0 = State.circular(Earth, a_0 * u.km - Earth.R, i_0 * u.rad)%0A    r0, v0 = s0.rv()%0A%0A    tf = 150 * s0.period%0A%0A    # Propagate orbit%0A    r, v = cowell(k,%0A                  r0.to(u.km).value,%0A                  v0.to(u.km / u.s).value,%0A                  tf.to(u.s).value,%0A                  ad=a_d,%0A                  nsteps=100000)%0A%0A    sf = State.from_vectors(Earth,%0A                            r * u.km,%0A                            v * u.km / u.s,%0A                            s0.epoch + tf)%0A%0A    print(sf.a.to(u.km))%0A    print(sf.ecc)%0A    print(sf.inc.to(%22deg%22))%0A

Add test to study the change in inclination only