Merge pull request #76 from hannorein/many_asteroids

Large asteroid file

MANIFEST.in

@@ -4,7 +4,6 @@ include src/planets.c
 include src/spk.c
 include src/assist.h
 include src/forces.h
-include src/const.h
 include src/planets.h
 include src/spk.h
 include README.md

docs/installation.md

@@ -1,5 +1,8 @@
 # Installation
 
+!!! Note 
+    If you have previously installed REBOUND in your environment, make sure you upgrade it to the latest version.
+
 === "Python"
     It's easiest to install ASSIST into a python virtual environment. If you already have a virtual environment or do not want to use one, you can skip this step. Otherwise, run the following command in an empty directory. They will setup and activate a new virtual environment in a directory. 
     ```bash
@@ -22,7 +25,7 @@
 
 !!! Warning inline end "Large file size"
 
-    Some of the files are over 2 GB in size.
+    Some of the files are over 2 GB in size. The sb441-n373.bsp file below is over 14 GB.
 
 To use ASSIST, you also need to download ephemeris data files. One file for planet ephemeris and another suplementary file for asteroid ephemeris. The following commands download these files with curl. You can also manually download them using your browser. You can store these files anywhere. If you're installing REBOUND from the repository, a good place to put them is in a new directory with the name `data`.
 
@@ -38,7 +41,10 @@ For some of the examples, you will also need the planet ephemeris file with an e
 curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441 -o assist/data/linux_m13000p17000.441
 ```
 
-If you wish to carry out a calculation with 373 small bodies (instead of just the 16 most massive asteroids in the  'n16' file above) you will want to
+!!! Info 
+    Support for the n373 file has been added recently. Compare your results to those obtained with the standard n16 file.
+
+If you wish to carry out a calculation with 373 small bodies (instead of just the 16 most massive asteroids in the  n16 file above) you will want to
 download the full asteroid file:
 
 ```bash
@@ -47,7 +53,6 @@ curl https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n373.bs
 
 You can now verify that your installation of assist works.
 
-Note: make sure you upgrade to the latest version of REBOUND if you have previously installed it in your envirnoment.
 
 === "Python"
     Start a python interpreter, or open a new notebook if you're planning to use Jupyter notebooks.

src/Makefile

@@ -24,7 +24,7 @@ endif
 
 SOURCES=assist.c spk.c planets.c forces.c
 OBJECTS=$(SOURCES:.c=.o)
-HEADERS=$(SOURCES:.c=.h) const.h
+HEADERS=$(SOURCES:.c=.h)
 
 all: $(SOURCES) librebound.so libassist.so
 

src/assist.c

@@ -23,8 +23,6 @@
  *
  */
 
-#include "const.h"
-
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
@@ -136,6 +134,81 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char
     if ((ephem->spl = assist_spk_init(asteroids_path)) == NULL) {
         return ASSIST_ERROR_AST_FILE;	  
     }
+    // Try to find masses of bodies in spk file in ephemeris constants
+    for(int n=0; n<ephem->spl->num; n++){ // loop over all asteroids
+        int found = 0;
+        for(int c=0; c<ephem->jpl->num; c++){ // loop over all constants
+	        if (strncmp(ephem->jpl->str[c], "MA", 2) == 0) {
+                int cid = atoi(ephem->jpl->str[c]+2);
+                int offset = 2000000;
+                if (cid==ephem->spl->targets[n].code-offset){
+                    ephem->spl->targets[n].mass = ephem->jpl->con[c];
+                    found = 1;
+                    break;
+                }
+            }
+        }
+        // Use lookup table for new KBO objects in DE440/441
+        // Source: https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/README.txt
+        int massmap[] = {
+            // ID, SPK_ID
+            8001,  2136199,
+            8002,  2136108,
+            8003,  2090377,
+            8004,  2136472,
+            8005,  2050000,
+            8006,  2084522,
+            8007,  2090482,
+            8008,  2020000,
+            8009,  2055637,
+            8010,  2028978,
+            8011,  2307261,
+            8012,  2174567,
+            8013,  3361580,
+            8014,  3308265,
+            8015,  2055565,
+            8016,  2145452,
+            8017,  2090568,
+            8018,  2208996,
+            8019,  2225088,
+            8020,  2019521,
+            8021,  2120347,
+            8022,  2278361,
+            8023,  3525142,
+            8024,  2230965,
+            8025,  2042301,
+            8026,  2455502,
+            8027,  3545742,
+            8028,  2523639,
+            8029,  2528381,
+            8030,  3515022,
+        };
+        if (found==0){
+            int mapped = -1;
+            for (int m=0; m<sizeof(massmap); m+=2){
+                if (massmap[m+1]==ephem->spl->targets[n].code){
+                    mapped = massmap[m];
+                    break;
+                }
+            }
+            if (mapped != -1){
+                for(int c=0; c<ephem->jpl->num; c++){ // loop over all constants (again)
+                    if (strncmp(ephem->jpl->str[c], "MA", 2) == 0) {
+                        int cid = atoi(ephem->jpl->str[c]+2);
+                        if (cid==mapped){
+                            ephem->spl->targets[n].mass = ephem->jpl->con[c];
+                            found = 1;
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+        if (found==0){
+            fprintf(stderr,"WARNING: Cannot find mass for asteroid %d (NAIF ID Number %d).\n", n, ephem->spl->targets[n].code );
+        }
+
+    }
 
     return ASSIST_SUCCESS;
 }
@@ -197,7 +270,7 @@ void assist_extras_cleanup(struct reb_simulation* sim){
 void assist_init(struct assist_extras* assist, struct reb_simulation* sim, struct assist_ephem* ephem){
     assist->sim = sim;
     assist->ephem_cache = calloc(1, sizeof(struct assist_ephem_cache));
-    const int N_total = ASSIST_BODY_NPLANETS + ASSIST_BODY_NASTEROIDS;
+    const int N_total = ASSIST_BODY_NPLANETS + ephem->spl->num;
     assist->gr_eih_sources = 1; // Only include Sun by default
     assist->ephem_cache->items = calloc(N_total*7, sizeof(struct assist_cache_item));
     assist->ephem_cache->t = malloc(N_total*7*sizeof(double));

src/assist.h

@@ -88,25 +88,6 @@ enum ASSIST_BODY {
     ASSIST_BODY_PLUTO       = 10,
 
     ASSIST_BODY_NPLANETS    = 11,
-
-    ASSIST_BODY_CAMILLA     = 11,
-    ASSIST_BODY_CERES       = 12,
-    ASSIST_BODY_CYBELE      = 13,
-    ASSIST_BODY_DAVIDA      = 14,
-    ASSIST_BODY_EUNOMIA     = 15,
-    ASSIST_BODY_EUPHROSYNE  = 16,
-    ASSIST_BODY_EUROPA      = 17,
-    ASSIST_BODY_HYGIEA      = 18,
-    ASSIST_BODY_INTERAMNIA  = 19,
-    ASSIST_BODY_IRIS        = 20,
-    ASSIST_BODY_JUNO        = 21,
-    ASSIST_BODY_PALLAS      = 22,
-    ASSIST_BODY_PSYCHE      = 23,
-    ASSIST_BODY_SYLVIA      = 24,
-    ASSIST_BODY_THISBE      = 25,
-    ASSIST_BODY_VESTA       = 26,
-
-    ASSIST_BODY_NASTEROIDS  = 16,
 };
 
 struct assist_ephem {