diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 1b140f119..181023676 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -63,9 +63,9 @@ def solar_system_horizons(name: str, ispl = False print(f"\nMassive body {name} not found or not yet supported") print("This will be created as a massless test particle") - if id is None: - print("ID value required for this input type") - return + # if id is None: + # print("ID value required for this input type") + # return # Planet MSun/M ratio MSun_over_Mpl = { @@ -177,9 +177,14 @@ def solar_system_horizons(name: str, ephemerides_step = '1d' pldata = {} - pldata[name] = Horizons(id=id, location='@sun', - epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, - 'step': ephemerides_step}) + if id is not None: + pldata[name] = Horizons(id=id, location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) + else: + pldata[name] = Horizons(id=name, location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) if param['IN_FORM'] == 'XV': rx = pldata[name].vectors()['x'][0] * DCONV @@ -221,7 +226,7 @@ def solar_system_horizons(name: str, else: Gmass = None - if id is None: + if id is None and name in planetid: id = planetid[name] return id,name,a,e,inc,capom,omega,capm,rh,vh,Gmass,Rpl,rhill,Ip,rot,J2,J4 diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 399ab0d38..f91a1871b 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2241,6 +2241,7 @@ def add_solar_system_body(self, for k,v in kwargs.items(): if k in scalar_ints: + v[v == None] = -1 kwargs[k] = v.astype(int) elif k in scalar_floats: kwargs[k] = v.astype(np.float64) @@ -2254,11 +2255,21 @@ def add_solar_system_body(self, kwargs['time'] = np.array([self.param['TSTART']]) + if len(self.data) == 0: + maxid = -1 + else: + maxid = self.data.id.max().values[()] + + nbodies = kwargs["name"].size + kwargs['id'] = np.where(kwargs['id'] < 0,np.arange(start=maxid+1,stop=maxid+1+nbodies,dtype=int),kwargs['id']) + dsnew = init_cond.vec2xr(self.param,**kwargs) dsnew = self._combine_and_fix_dsnew(dsnew) if dsnew['npl'] > 0 or dsnew['ntp'] > 0: self.save(verbose=False) + + self.init_cond = self.data.copy(deep=True) diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index fdcd5585f..e81194d37 100644 --- a/python/swiftest/swiftest/tool.py +++ b/python/swiftest/swiftest/tool.py @@ -359,7 +359,7 @@ def el2xv_vec(mu, a, ecc, inc, Omega, omega, M): longitude of ascending node (degrees) omega : (n) float array argument of periapsis (degrees) - M : (n) flaot array + M : (n) float array Mean anomaly (degrees) Returns