From 848cecc10babc8b0f8d40c42a4fab4c2a551aa87 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 23:01:19 -0500 Subject: [PATCH] Added GR test script and Notebook to helio_gr_test --- .../helio_gr_test/swiftest_relativity.ipynb | 177 ++++++++++++++++++ examples/helio_gr_test/swiftest_relativity.py | 60 ++++++ 2 files changed, 237 insertions(+) create mode 100644 examples/helio_gr_test/swiftest_relativity.ipynb create mode 100644 examples/helio_gr_test/swiftest_relativity.py diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb new file mode 100644 index 000000000..6abb9524a --- /dev/null +++ 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=\"whm\",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=\"whm\",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 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", + "plt.savefig(\"whm_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 new file mode 100644 index 000000000..a5f4e4371 --- /dev/null +++ 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)}')