Altered example for J2 test case so that it can compare the c_lm vs J…

…2 versions. Also fixed some issues with getting central body initial conditions in correctly when in EL input mode

examples/spherical_harmonics_cb/J2_test_tp.py

@@ -24,9 +24,9 @@
 seed = 123
 rng = np.random.default_rng(seed=seed)
 
-# set up swiftest simulation with relevant units (here they are km, days, and kg)
-sim = swiftest.Simulation(DU2M = 1e3, TU = 'd', MU = 'kg')
-sim.clean()
+
+
+
 
 # Central Body Parameters (Chariklo parameters from Leiva, et al (2017) (Jacobi Ellipsoid model))
 cb_mass = 6.1e18 # kg
@@ -39,6 +39,17 @@
 cb_T_rotation = 7.004 / 24.0 # converting from hours to julian days (TU)
 cb_rot = [[0, 0, 360.0 / cb_T_rotation]] # degrees/d
 
+# Add 1 user-defined test particle.
+ntp = 1
+
+name_tp     = ["TestParticle_01"]
+a_tp        = 300
+e_tp        = 0.05
+inc_tp      = 10
+capom_tp    = 0.0
+omega_tp    = 0.0
+capm_tp     = 0.0
+
 # Extract the spherical harmonics coefficients (c_lm) from axes measurements
 #
 # The user can pass an optional reference radius at which the coefficients are calculated. If not provided, SHTOOLS 
@@ -55,26 +66,21 @@
 J2 = -tmp20 * np.sqrt(5) # unnormalised J2 term
 j2rp2 = J2 * cb_radius**2
 
-# Add the central body
-# The user can pass the c_lm coefficients directly to the add_body method if they do not wish to use the clm_from_ellipsoid method.
-sim.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, c_lm = c_lm)
-
-# Add 1 user-defined test particle.
-ntp = 1
+# set up swiftest simulation with relevant units (here they are km, days, and kg)
+sim_shgrav = swiftest.Simulation(simdir="shgrav",DU2M = 1e3, TU = 'd', MU = 'kg')
 
-name_tp     = ["TestParticle_01"]
-a_tp        = 300
-e_tp        = 0.05
-inc_tp      = 10
-capom_tp    = 0.0
-omega_tp    = 0.0
-capm_tp     = 0.0
+sim_shgrav.clean()
+# Use the shgrav version where you input a set of spherical harmonics coefficients
+sim_shgrav.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, c_lm = c_lm)
+sim_shgrav.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
+sim_shgrav.run(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True)
 
-sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
-sim.set_parameter(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True)
+# Use the original "oblate" version where you pass J2 (and/or J4)
+sim_obl = swiftest.Simulation(simdir="obl", DU2M = 1e3, TU='d', MU='kg')
+sim_obl.clean()
+sim_obl.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, J2 = j2rp2)
+sim_obl.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
+sim_obl.run(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True)
 
-# Display the run configuration parameters.
-sim.get_parameter()
+ds_diff = sim_shgrav.data - sim_obl.data
 
-# Run the simulation. Arguments may be defined here or thorugh the swiftest.Simulation() method.
-sim.run()

swiftest/simulation_class.py

@@ -2637,17 +2637,18 @@ def input_to_clm_array(val, n):
                 raise ValueError("Orbital elements cannot be passed for a central body.")
             if nbodies > 1:
                 raise ValueError("Only one central body may be passed.")
-            if rh is None:
-                rh = np.zeros((1,3))
-            if vh is None:
-                vh = np.zeros((1,3))
-            a = np.nan
-            e = np.nan
-            inc = np.nan
-            capom = np.nan
-            omega = np.nan
-            capm = np.nan
-
+            if self.param['IN_FORM'] == "XV":
+                if rh is None:
+                    rh = np.zeros((1,3))
+                if vh is None:
+                    vh = np.zeros((1,3))
+            elif self.param['IN_FORM'] == "EL":
+                a = np.array([np.nan])
+                e = np.array([np.nan])
+                inc = np.array([np.nan])
+                capom = np.array([np.nan])
+                omega = np.array([np.nan])
+                capm = np.array([np.nan])
 
         dsnew = init_cond.vec2xr(self.param, name=name, a=a, e=e, inc=inc, capom=capom, omega=omega, capm=capm, id=id,
                                  Gmass=Gmass, radius=radius, rhill=rhill, Ip=Ip, rh=rh, vh=vh,rot=rot, j2rp2=J2, j4rp4=J4, c_lm=c_lm, rotphase=rotphase, time=time)