updated whm_gr_test

MintonGroup · Dec 4, 2022 · eadc6a6 · eadc6a6
1 parent d972a21
commit eadc6a6
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 41 deletions.
diff --git a/...les/whm_gr_test/swiftest_relativity.ipynb → examples/whm_gr_test/whm_gr_test.ipynb b/...les/whm_gr_test/swiftest_relativity.ipynb → examples/whm_gr_test/whm_gr_test.ipynb
@@ -10,8 +10,7 @@
     "from astroquery.jplhorizons import Horizons\n",
     "import datetime\n",
     "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "import os"
+    "import matplotlib.pyplot as plt"
    ]
   },
   {
@@ -89,17 +88,6 @@
     "varpi_obs = el['w'] + el['Omega']"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Compute the longitude of the periapsis\n",
-    "sim_gr.data['varpi'] = np.mod(sim_gr.data['omega'] + sim_gr.data['capom'],360)\n",
-    "sim_nogr.data['varpi'] = np.mod(sim_nogr.data['omega'] + sim_nogr.data['capom'],360)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -131,8 +119,8 @@
     "fig, ax = plt.subplots()\n",
     "\n",
     "ax.plot(t, varpi_obs, label=\"JPL Horizons\",linewidth=2.5)\n",
-    "ax.plot(tsim, varpisim_gr, label=\"Swiftest whm GR\",linewidth=1.5)\n",
-    "ax.plot(tsim, varpisim_nogr, label=\"Swiftest whm No GR\",linewidth=1.5)\n",
+    "ax.plot(tsim, varpisim_gr, label=\"Swiftest WHM GR\",linewidth=1.5)\n",
+    "ax.plot(tsim, varpisim_nogr, label=\"Swiftest WHM No GR\",linewidth=1.5)\n",
     "ax.set_xlabel('Time (y)')\n",
     "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n",
     "ax.legend()\n",
@@ -144,22 +132,6 @@
     "print(f'Obs - Swiftest GR    : {np.mean(dvarpi_obs - dvarpi_gr)}')\n",
     "print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sim_nogr.param"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

diff --git a/examples/whm_gr_test/swiftest_relativity.py → examples/whm_gr_test/whm_gr_test.py b/examples/whm_gr_test/swiftest_relativity.py → examples/whm_gr_test/whm_gr_test.py
@@ -11,9 +11,10 @@
 sim_nogr = swiftest.Simulation(simdir="nogr")
 sim_nogr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
 
-tstep_out = 10.0
-sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="whm",general_relativity=True)
-sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="whm",general_relativity=False)
+run_args = {"tstop":1000.0, "dt":0.005, "tstep_out":10.0, "dump_cadence": 0,"integrator":"whm"}
+
+sim_gr.run(**run_args,general_relativity=True)
+sim_nogr.run(**run_args,general_relativity=False)
 
 # Get the start and end date of the simulation so we can compare with the real solar system
 start_date = sim_gr.ephemeris_date
@@ -29,19 +30,14 @@
 t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25
 varpi_obs = el['w'] + el['Omega']
 
-# Compute the longitude of the periapsis
-sim_gr.data['varpi'] = np.mod(sim_gr.data['omega'] + sim_gr.data['capom'],360)
-sim_nogr.data['varpi'] = np.mod(sim_nogr.data['omega'] + sim_nogr.data['capom'],360)
-
 varpisim_gr= sim_gr.data['varpi'].sel(name="Mercury")
 varpisim_nogr= sim_nogr.data['varpi'].sel(name="Mercury")
 tsim = sim_gr.data['time']
 
-dvarpi_gr = np.diff(varpisim_gr)  * 3600 * 100 / tstep_out
-dvarpi_nogr = np.diff(varpisim_nogr)  * 3600 * 100 / tstep_out
+dvarpi_gr = np.diff(varpisim_gr)  * 3600 * 100 / run_args['tstep_out']
+dvarpi_nogr = np.diff(varpisim_nogr)  * 3600 * 100 / run_args['tstep_out']
 dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100
 
-
 fig, ax = plt.subplots()
 
 ax.plot(t, varpi_obs, label="JPL Horizons",linewidth=2.5)