diff --git a/examples/helio_gr_test/grsim/param.gr.in b/examples/helio_gr_test/grsim/param.gr.in deleted file mode 100644 index 0616db203..000000000 --- a/examples/helio_gr_test/grsim/param.gr.in +++ /dev/null @@ -1,35 +0,0 @@ -! VERSION Swiftest input file -T0 0.0 -TSTART 0.0 -TSTOP 1000.0 -DT 0.005 -ISTEP_OUT 2000 -ISTEP_DUMP 2000 -NC_IN init_cond.nc -IN_TYPE NETCDF_DOUBLE -IN_FORM EL -BIN_OUT bin.gr.nc -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -CHK_QMIN 0.004650467260962157 -CHK_RMIN 0.004650467260962157 -CHK_RMAX 10000.0 -CHK_EJECT 10000.0 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.004650467260962157 10000.0 -MU2KG 1.988409870698051e+30 -TU2S 31557600.0 -DU2M 149597870700.0 -FRAGMENTATION NO -RESTART NO -CHK_CLOSE YES -GR YES -ROTATION NO -ENERGY NO -EXTRA_FORCE NO -BIG_DISCARD NO -RHILL_PRESENT NO -INTERACTION_LOOPS TRIANGULAR -ENCOUNTER_CHECK TRIANGULAR -TIDES NO diff --git a/examples/helio_gr_test/nogrsim/param.nogr.in b/examples/helio_gr_test/nogrsim/param.nogr.in deleted file mode 100644 index 9e2ab0b22..000000000 --- a/examples/helio_gr_test/nogrsim/param.nogr.in +++ /dev/null @@ -1,35 +0,0 @@ -! VERSION Swiftest input file -T0 0.0 -TSTART 0.0 -TSTOP 1000.0 -DT 0.005 -ISTEP_OUT 2000 -ISTEP_DUMP 2000 -NC_IN init_cond.nc -IN_TYPE NETCDF_DOUBLE -IN_FORM EL -BIN_OUT bin.nogr.nc -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -CHK_QMIN 0.004650467260962157 -CHK_RMIN 0.004650467260962157 -CHK_RMAX 10000.0 -CHK_EJECT 10000.0 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.004650467260962157 10000.0 -MU2KG 1.988409870698051e+30 -TU2S 31557600.0 -DU2M 149597870700.0 -FRAGMENTATION NO -RESTART NO -CHK_CLOSE YES -GR NO -ROTATION NO -ENERGY NO -EXTRA_FORCE NO -BIG_DISCARD NO -RHILL_PRESENT NO -INTERACTION_LOOPS TRIANGULAR -ENCOUNTER_CHECK TRIANGULAR -TIDES NO diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb deleted file mode 100644 index e7ae73b23..000000000 --- a/examples/helio_gr_test/swiftest_relativity.ipynb +++ /dev/null @@ -1,176 +0,0 @@ -{ - "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(simdir=\"gr\")\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(simdir=\"nogr\")\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": [ - "run_args = {\"tstop\":1000.0, \"dt\":0.005, \"tstep_out\":10.0, \"dump_cadence\": 0,\"integrator\":\"helio\"}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sim_gr.run(**run_args,general_relativity=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sim_nogr.run(**run_args,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 / run_args['tstep_out']\n", - "dvarpi_nogr = np.diff(varpisim_nogr) * 3600 * 100 / run_args['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": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python (My debug_env Kernel)", - "language": "python", - "name": "debug_env" - }, - "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 deleted file mode 100644 index a5f4e4371..000000000 --- a/examples/helio_gr_test/swiftest_relativity.py +++ /dev/null @@ -1,60 +0,0 @@ -#!/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)}') diff --git a/examples/whm_gr_test/swiftest_relativity.py b/examples/whm_gr_test/swiftest_relativity.py index a4ea53c3b..d4b777fa1 100644 --- a/examples/whm_gr_test/swiftest_relativity.py +++ b/examples/whm_gr_test/swiftest_relativity.py @@ -5,10 +5,10 @@ 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 = swiftest.Simulation(simdir="gr") 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 = swiftest.Simulation(simdir="nogr") sim_nogr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"]) tstep_out = 10.0