From f6e1a9e5e1229d05170dc81574df200b5ed97129 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 14 Dec 2021 10:59:02 -0500 Subject: [PATCH] Fixed bad array index range value (should be ntp not npl) --- examples/symba_gr_test/init_cond.py | 27 ++++++++++++++++++++++++++- src/symba/symba_step.f90 | 2 +- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/examples/symba_gr_test/init_cond.py b/examples/symba_gr_test/init_cond.py index 2639272e7..a5dbe9118 100755 --- a/examples/symba_gr_test/init_cond.py +++ b/examples/symba_gr_test/init_cond.py @@ -1,5 +1,7 @@ #!/usr/bin/env python3 import swiftest +from numpy.random import default_rng +import numpy as np sim = swiftest.Simulation() sim.param['PL_IN'] = "pl.swiftest.in" @@ -43,7 +45,30 @@ for name, id in bodyid.items(): sim.add(name, idval=id, date="2027-04-30") - + +Me_a = sim.ds.isel(id=1)['a'].values +Me_e = sim.ds.sel(id=1)['e'].values +Me_i = sim.ds.sel(id=1)['inc'].values + +capom_pl = default_rng().uniform(0.0, 360.0, 1) +omega_pl = default_rng().uniform(0.0, 360.0, 1) +capm_pl = default_rng().uniform(0.0, 360.0, 1) + +capom_tp = default_rng().uniform(0.0, 360.0, 1) +omega_tp = default_rng().uniform(0.0, 360.0, 1) +capm_tp = default_rng().uniform(0.0, 360.0, 1) + +GMcb = sim.ds.isel(id=0)['Gmass'].values +GU = swiftest.GC / (sim.param['DU2M']**3 / (sim.param['MU2KG'] * sim.param['TU2S']**2)) +dens = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M']**3) # Assume a bulk density of 3 g/cm^3 +GM_pl = 2e-7 +M_pl = GM_pl / GU +R_pl = (3 * M_pl / (4 * np.pi * dens))**(1.0 / 3.0) +Rh_pl = Me_a * (GM_pl / (3 * GMcb))**(1.0/3.0) + +sim.addp(np.full(1,9), np.full(1,'Planetesimal'), Me_a, Me_e, Me_i, capom_pl, omega_pl, capm_pl, GMpl=np.full(1, GM_pl), Rpl=np.full(1, R_pl), rhill=Rh_pl) +sim.addp(np.full(1,10), np.full(1,'TestParticle'), Me_a, Me_e, Me_i, capom_tp, omega_tp, capm_tp) + sim.save("param.swiftest.in") sim.param['PL_IN'] = "pl.swifter.in" sim.param['TP_IN'] = "tp.swifter.in" diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 2bd6aca36..cd48667af 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -274,7 +274,7 @@ module subroutine symba_step_reset_system(self, param) tp%levelg(1:ntp) = -1 tp%levelm(1:ntp) = -1 tp%lmask(1:ntp) = .true. - tp%ldiscard(1:npl) = .false. + tp%ldiscard(1:ntp) = .false. nenc_old = system%pltpenc_list%nenc call system%pltpenc_list%setup(0) call system%pltpenc_list%setup(nenc_old)