diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb new file mode 100644 index 000000000..6946ef658 --- /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=\"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 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)}') diff --git a/examples/whm_gr_test/init_cond.py b/examples/whm_gr_test/init_cond.py deleted file mode 100644 index 7904eb100..000000000 --- a/examples/whm_gr_test/init_cond.py +++ /dev/null @@ -1,226 +0,0 @@ -import numpy as np -import sys -from astroquery.jplhorizons import Horizons -import astropy.constants as const - -#Values from JPL Horizons -AU2M = const.au.value -GMSunSI = const.GM_sun.value -Rsun = const.R_sun.value -GC = const.G.value -JD = 86400 -year = 365.25 * JD -c = 299792458.0 -MSun_over_Mpl = [6023600.0, - 408523.71, - 328900.56, - 3098708., - 1047.3486, - 3497.898, - 22902.98, - 19412.24, - 1.35e8] - -MU2KG = GMSunSI / GC #Conversion from mass unit to kg -DU2M = AU2M #Conversion from radius unit to centimeters -TU2S = year #Conversion from time unit to seconds -GU = GC / (DU2M**3 / (MU2KG * TU2S**2)) - -GMSun = GMSunSI / (DU2M**3 / TU2S**2) - -t_print = 10.e0 * year / TU2S #output interval to print results -deltaT = 0.25 * JD / TU2S #timestep simulation -end_sim = 1.0e3 * year / TU2S + t_print #end time - -# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) -J2 = 2.198e-7 * (Rsun / DU2M)**2 -J4 = -4.805e-9 * (Rsun / DU2M)**4 - -tstart = '2021-01-28' -tend = '2021-01-29' -tstep = '1d' -planetid = { - 'mercury' : '1', - 'venus' : '2', - 'earthmoon' : '3', - 'mars' : '4', - 'jupiter' : '5', - 'saturn' : '6', - 'uranus' : '7', - 'neptune' : '8', - 'plutocharon' : '9' -} -npl = 9 - -#Planet Msun/M ratio -MSun_over_Mpl = { - 'mercury' : 6023600.0, - 'venus' : 408523.71, - 'earthmoon' : 328900.56, - 'mars' : 3098708., - 'jupiter' : 1047.3486, - 'saturn' : 3497.898, - 'uranus' : 22902.98, - 'neptune' : 19412.24, - 'plutocharon' : 1.35e8 -} - -#Planet radii in meters -Rpl = { - 'mercury' : 2439.4e3, - 'venus' : 6051.8e3, - 'earthmoon' : 6371.0084e3, # Earth only for radius - 'mars' : 3389.50e3, - 'jupiter' : 69911e3, - 'saturn' : 58232.0e3, - 'uranus' : 25362.e3, - 'neptune' : 24622.e3, - 'plutocharon' : 1188.3e3 -} - -pdata = {} -plvec = {} -Rhill = {} - -for key,val in planetid.items(): - pdata[key] = Horizons(id=val, id_type='majorbody',location='@sun', - epochs={'start': tstart, 'stop': tend, - 'step': tstep}) - plvec[key] = np.array([pdata[key].vectors()['x'][0], - pdata[key].vectors()['y'][0], - pdata[key].vectors()['z'][0], - pdata[key].vectors()['vx'][0], - pdata[key].vectors()['vy'][0], - pdata[key].vectors()['vz'][0] - ]) - Rhill[key] = pdata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0) - - -if __name__ == '__main__': - # Convert from AU-day to AU-year just because I find it easier to keep track of the sim progress - for plid in plvec: - plvec[plid][3:] *= year / JD - - # Names of all output files - swifter_input = "param.swifter.in" - swifter_pl = "pl.swifter.in" - swifter_tp = "tp.swifter.in" - swifter_bin = "bin.swifter.dat" - swifter_enc = "enc.swifter.dat" - - swiftest_input = "param.swiftest.in" - swiftest_pl = "pl.swiftest.in" - swiftest_tp = "tp.swiftest.in" - swiftest_cb = "cb.swiftest.in" - swiftest_bin = "bin.swiftest.dat" - swiftest_enc = "enc.swiftest.dat" - - # Simulation start, stop, and output cadence times - t_0 = 0 # simulation start time - end_sim = 1000.0e0 * year / TU2S # simulation end time - deltaT = 0.25 * JD / TU2S # simulation step size - t_print = 1.0 * year / TU2S #output interval to print results - - iout = int(np.ceil(t_print / deltaT)) - rmin = Rsun / DU2M - rmax = 1000.0 - - #Make Swifter files - plfile = open(swifter_pl, 'w') - print(f'{npl+1} ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch {tstart}' ,file=plfile) - print(f'1 {GMSun}',file=plfile) - print(f'0.0 0.0 0.0',file=plfile) - print(f'0.0 0.0 0.0',file=plfile) - for i, plid in enumerate(plvec): - print(f'{i + 2} {GMSun / MSun_over_Mpl[plid]} {Rhill[plid]}', file=plfile) - print(f'{Rpl[plid] / DU2M}', file=plfile) - print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) - print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) - plfile.close() - - tpfile = open(swifter_tp, 'w') - print(0,file=tpfile) - tpfile.close() - - sys.stdout = open(swifter_input, "w") - print(f'! Swifter input file generated using init_cond.py') - print(f'T0 {t_0} ') - print(f'TSTOP {end_sim}') - print(f'DT {deltaT}') - print(f'PL_IN {swifter_pl}') - print(f'TP_IN {swifter_tp}') - print(f'IN_TYPE ASCII') - print(f'ISTEP_OUT {iout:d}') - print(f'ISTEP_DUMP {iout:d}') - print(f'BIN_OUT {swifter_bin}') - print(f'OUT_TYPE REAL8') - print(f'OUT_FORM EL') - print(f'OUT_STAT NEW') - print(f'J2 {J2}') - print(f'J4 {J4}') - print(f'CHK_CLOSE yes') - print(f'CHK_RMIN {rmin}') - print(f'CHK_RMAX {rmax}') - print(f'CHK_EJECT {rmax}') - print(f'CHK_QMIN {rmin}') - print(f'CHK_QMIN_COORD HELIO') - print(f'CHK_QMIN_RANGE {rmin} {rmax}') - print(f'ENC_OUT {swifter_enc}') - print(f'EXTRA_FORCE no') - print(f'BIG_DISCARD no') - print(f'RHILL_PRESENT yes') - print(f'C {c / (DU2M / TU2S)}') - - #Now make Swiftest files - cbfile = open(swiftest_cb, 'w') - print(f'{1.0}',file=cbfile) - print(f'{rmin}',file=cbfile) - print(f'{J2}',file=cbfile) - print(f'{J4}',file=cbfile) - - plfile = open(swiftest_pl, 'w') - print(npl,file=plfile) - - for i, plid in enumerate(plvec): - print(f'{i + 2} {1.0 / MSun_over_Mpl[plid]}', file=plfile) - print(f'{Rpl[plid] / DU2M}', file=plfile) - print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) - print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) - plfile.close() - tpfile = open(swiftest_tp, 'w') - print(0,file=tpfile) - tpfile.close() - - sys.stdout = open(swiftest_input, "w") - print(f'! Swiftest input file generated using init_cond.py') - print(f'T0 {t_0} ') - print(f'TSTOP {end_sim}') - print(f'DT {deltaT}') - print(f'CB_IN {swiftest_cb}') - print(f'PL_IN {swiftest_pl}') - print(f'TP_IN {swiftest_tp}') - print(f'IN_TYPE ASCII') - print(f'ISTEP_OUT {iout:d}') - print(f'ISTEP_DUMP {iout:d}') - print(f'BIN_OUT {swiftest_bin}') - print(f'OUT_TYPE REAL8') - print(f'OUT_FORM EL') - print(f'OUT_STAT REPLACE') - print(f'CHK_CLOSE yes') - print(f'CHK_RMIN {rmin}') - print(f'CHK_RMAX {rmax}') - print(f'CHK_EJECT {rmax}') - print(f'CHK_QMIN {rmin}') - print(f'CHK_QMIN_COORD HELIO') - print(f'CHK_QMIN_RANGE {rmin} {rmax}') - print(f'ENC_OUT {swiftest_enc}') - print(f'EXTRA_FORCE no') - print(f'BIG_DISCARD no') - print(f'ROTATION no') - print(f'GR yes') - print(f'MU2KG {MU2KG}') - print(f'DU2M {DU2M}') - print(f'TU2S {TU2S}') - - - sys.stdout = sys.__stdout__ diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index 6abb9524a..113e10f81 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -15,9 +15,448 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Sun as a central body\n", + "Fetching ephemerides data for Mercury from JPL/Horizons\n", + "Fetching ephemerides data for Venus from JPL/Horizons\n", + "Fetching ephemerides data for Earth from JPL/Horizons\n", + "Fetching ephemerides data for Mars from JPL/Horizons\n", + "Fetching ephemerides data for Jupiter from JPL/Horizons\n", + "Fetching ephemerides data for Saturn from JPL/Horizons\n", + "Fetching ephemerides data for Uranus from JPL/Horizons\n", + "Fetching ephemerides data for Neptune from JPL/Horizons\n", + "Writing initial conditions to file init_cond.nc\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/whm_gr_test/param.gr.in\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    a              (time, name) float64 nan 0.3871 0.7233 ... 9.532 19.24 30.04\n",
+       "    e              (time, name) float64 nan 0.2056 0.006718 ... 0.04796 0.008956\n",
+       "    inc            (time, name) float64 nan 7.003 3.394 ... 2.488 0.773 1.771\n",
+       "    capom          (time, name) float64 nan 48.3 76.6 ... 113.6 74.01 131.8\n",
+       "    ...             ...\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    a              (time, name) float64 nan 0.3871 0.7233 ... 9.532 19.24 30.04\n",
+       "    e              (time, name) float64 nan 0.2056 0.006718 ... 0.04796 0.008956\n",
+       "    inc            (time, name) float64 nan 7.003 3.394 ... 2.488 0.773 1.771\n",
+       "    capom          (time, name) float64 nan 48.3 76.6 ... 113.6 74.01 131.8\n",
+       "    ...             ...\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name)