Merge branch 'debug'

MintonGroup · Nov 18, 2022 · c8be13b · c8be13b
2 parents 0627031 + 60cf506
commit c8be13b
Show file tree

Hide file tree

Showing 5 changed files with 1,183 additions and 234 deletions.
diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb
@@ -0,0 +1,177 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import swiftest\n",
+    "from astroquery.jplhorizons import Horizons\n",
+    "import datetime\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim_gr = swiftest.Simulation(param_file=\"param.gr.in\", output_file_name=\"bin.gr.nc\")\n",
+    "sim_gr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim_nogr = swiftest.Simulation(param_file=\"param.nogr.in\", output_file_name=\"bin.nogr.nc\")\n",
+    "sim_nogr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "tstep_out = 10.0\n",
+    "sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%capture\n",
+    "sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get the start and end date of the simulation so we can compare with the real solar system\n",
+    "start_date = sim_gr.ephemeris_date\n",
+    "tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S\n",
+    "\n",
+    "stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Get the ephemerides of Mercury for the same timeframe as the simulation\n",
+    "obj = Horizons(id='1', location='@sun',\n",
+    "               epochs={'start':start_date, 'stop':stop_date,\n",
+    "                       'step':'10y'})\n",
+    "el = obj.elements()\n",
+    "t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25\n",
+    "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,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "varpisim_gr= sim_gr.data['varpi'].sel(name=\"Mercury\")\n",
+    "varpisim_nogr= sim_nogr.data['varpi'].sel(name=\"Mercury\")\n",
+    "tsim = sim_gr.data['time']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dvarpi_gr = np.diff(varpisim_gr)  * 3600 * 100 / tstep_out\n",
+    "dvarpi_nogr = np.diff(varpisim_nogr)  * 3600 * 100 / tstep_out\n",
+    "dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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 helio GR\",linewidth=1.5)\n",
+    "ax.plot(tsim, varpisim_nogr, label=\"Swiftest helio No GR\",linewidth=1.5)\n",
+    "ax.set_xlabel('Time (y)')\n",
+    "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n",
+    "ax.legend()\n",
+    "plt.savefig(\"helio_gr_mercury_precession.png\",dpi=300)\n",
+    "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n",
+    "print(f'JPL Horizons         : {np.mean(dvarpi_obs)}')\n",
+    "print(f'Swiftest No GR       : {np.mean(dvarpi_nogr)}')\n",
+    "print(f'Swiftest GR          : {np.mean(dvarpi_gr)}')\n",
+    "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": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "swiftest",
+   "language": "python",
+   "name": "swiftest"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/helio_gr_test/swiftest_relativity.py b/examples/helio_gr_test/swiftest_relativity.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+import swiftest
+from astroquery.jplhorizons import Horizons
+import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+
+sim_gr = swiftest.Simulation(param_file="param.gr.in", output_file_name="bin.gr.nc")
+sim_gr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
+
+sim_nogr = swiftest.Simulation(param_file="param.nogr.in", output_file_name="bin.nogr.nc")
+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="helio",general_relativity=True)
+sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="helio",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
+tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S
+
+stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat()
+
+#Get the ephemerides of Mercury for the same timeframe as the simulation
+obj = Horizons(id='1', location='@sun',
+               epochs={'start':start_date, 'stop':stop_date,
+                       'step':'10y'})
+el = obj.elements()
+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_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)
+ax.plot(tsim, varpisim_gr, label="Swiftest Helio GR",linewidth=1.5)
+ax.plot(tsim, varpisim_nogr, label="Swiftest Helio No GR",linewidth=1.5)
+ax.set_xlabel('Time (y)')
+ax.set_ylabel('Mercury $\\varpi$ (deg)')
+ax.legend()
+plt.savefig("helio_gr_mercury_precession.png",dpi=300)
+
+print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')
+print(f'JPL Horizons         : {np.mean(dvarpi_obs)}')
+print(f'Swiftest No GR       : {np.mean(dvarpi_nogr)}')
+print(f'Swiftest GR          : {np.mean(dvarpi_gr)}')
+print(f'Obs - Swiftest GR    : {np.mean(dvarpi_obs - dvarpi_gr)}')
+print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')