Fixed add_solar_system_body method to be able to handle small bodies …

…properly

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

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)
 

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