examples updates

examples/test_GTO2GEO_caseB.py

@@ -0,0 +1,113 @@
+"""
+test with case B from Varga and Perez, 2016.
+"""
+
+
+import numpy as np
+from numpy.random import rand
+import matplotlib.pyplot as plt
+import time
+import spiceypy as spice
+
+import sys
+sys.path.append("../")
+import pyqlaw
+
+import os
+spice.furnsh(os.path.join(os.getenv("SPICE"), "lsk", "naif0012.tls"))
+spice.furnsh(os.path.join(os.getenv("SPICE"), "spk", "de440.bsp"))
+spice.furnsh(os.path.join(os.getenv("SPICE"), "pck", "gm_de440.tpc"))
+spice.furnsh(os.path.join(os.getenv("SPICE"), "pck", "earth_200101_990825_predict.bpc"))
+spice.furnsh(os.path.join(os.getenv("SPICE"), "fk", "earth_assoc_itrf93.tf"))
+
+plt.rcParams.update({'font.size': 15})
+
+def test_object():
+    tstart = time.time()
+    # reference quantities
+    MU_EARTH = spice.bodvrd("399", "GM", 1)[1][0]
+    LU = 42165.0
+    VU = np.sqrt(MU_EARTH / LU)
+    TU = LU / VU
+
+    # initialize object for perturbations
+    et_ref = spice.str2et("2028-01-01T00:00:00")
+    perturbations = pyqlaw.SpicePerturbations(et_ref, LU, TU)
+
+    # construct problem
+    prob = pyqlaw.QLaw(
+        integrator="rk4", 
+        elements_type="mee_with_a",
+        verbosity=2,
+        print_frequency=3000,
+        use_sundman = True,
+        perturbations = perturbations,
+        relaxed_tol_factor = 1,
+    )
+    # initial and final elements: [a,e,i,RAAN,omega,ta]
+    KEP0 = np.array([
+        24505.9/LU,
+        0.725,
+        np.deg2rad(7.05),
+        0.0, np.deg2rad(270), 0.0
+    ])
+    KEPF = np.array([
+        42165.0/LU,
+        0.001,
+        np.deg2rad(0.05),
+        np.deg2rad(10), np.deg2rad(10), np.deg2rad(10),
+    ])
+    oe0 = pyqlaw.kep2mee_with_a(np.array(KEP0))
+    oeT = pyqlaw.kep2mee_with_a(np.array(KEPF))
+    print(f"oe0: {oe0}")
+    print(f"oeT: {oeT}")
+    woe = [1.0, 1.0, 1.0, 1.0, 1.0]
+
+    # spacecraft parameters
+    MU = 2000
+    tmax_si = 1.0 #0.35   # N
+    isp_si  = 2000   # seconds
+    mdot_si = tmax_si/(isp_si*9.81)  # kg/s
+
+    # non-dimensional quantities
+    mass0 = 1.0
+    tmax = tmax_si * (1/MU)*(TU**2/(1e3*LU))
+    mdot = np.abs(mdot_si) *(TU/MU)
+    tf_max = 10000.0
+    t_step = np.deg2rad(2)
+
+    # set problem
+    prob.set_problem(oe0, oeT, mass0, tmax, mdot, tf_max, t_step, woe=woe)
+    prob.pretty()
+
+    # solve
+    prob.solve(eta_a=0.2, eta_r=0.3)
+    prob.pretty_results()
+    tend = time.time()
+    print(f"Simulation took {tend-tstart:4.4f} seconds")
+
+    # plot
+    fig1, ax1 = prob.plot_elements_history(to_keplerian=False)
+    fig2, ax2 = prob.plot_trajectory_3d(sphere_radius=0.1)
+    ax2.set(xlabel="x, LU", ylabel="y, LU", zlabel="z, LU")
+    fig3, ax3 = prob.plot_controls()
+    fig4, ax4 = prob.plot_efficiency()
+    fig5, ax5 = prob.plot_Q()
+
+    # # save
+    # fig2.savefig(os.path.join(os.path.dirname(__file__),
+    #                          "../paper/example_3D_trajectory.png"))
+    # fig1.savefig(os.path.join(os.path.dirname(__file__),
+    #                          "../paper/example_state_history.png"), bbox_inches='tight')
+    # fig3.savefig(os.path.join(os.path.dirname(__file__),
+    #                          "../paper/example_control_history.png"), bbox_inches='tight')
+
+    # export state history as initial guess for ICLOCS2
+    #prob.save_to_dict('initial_guess_GTO2GEO.json')
+    return fig1, fig2, fig3
+
+
+if __name__=="__main__":
+    figs = test_object()
+    plt.show()
+    print("Done!")

examples/test_GTO2GEO_perturbed.py

@@ -65,7 +65,7 @@ def test_object():
 
     # spacecraft parameters
     mass0 = 1.0
-    tmax = 0.0149
+    tmax = 0.0149/2
     mdot = 0.0031
     tf_max = 10000.0
     t_step = np.deg2rad(2)

examples/test_solve_mee_interplanetary.py

@@ -38,7 +38,7 @@ def test_object():
     tmax = 1e-2
     mdot = 1e-3
     tf_max = 2000.0
-    t_step = 0.1
+    t_step = 0.01
     # set problem
     prob.set_problem(oe0, oeT, mass0, tmax, mdot, tf_max, t_step, woe=woe)
     prob.pretty()

pyqlaw/_qlaw.py

@@ -666,7 +666,7 @@ def get_cartesian_history(self, interpolate=False, steps=None, kind="quadratic",
             # interpolate orbital elements
             f_a, f_e, f_i, f_r, f_o, f_t = self.interpolate_states(kind=kind)
             if steps is None:
-                steps = min(8000, abs(int(round(self.times[-1]/0.02))))
+                steps = min(8000, abs(int(round(self.times[-1]/0.01))))
                 print(f"Using {steps} steps for evaluation")
             t_evals = np.linspace(self.times[0], self.times[-1], steps)
             cart = np.zeros((6,steps))