Merge pull request #329 from VirtualPlanetaryLaboratory/StopAge

Stop age

Author: RoryBarnes

src/evolve.h

@@ -14,9 +14,13 @@
 #define EULER 1
 #define RUNGEKUTTA 2
 
+#define NO_INTEGRATION 0
+#define FORWARD_INTEGRATION 1
+#define BACKWARD_INTEGRATION -1
+
 /* @cond DOXYGEN_OVERRIDE */
 
-double fdTrapezoidalArea(double,double,double);
+double fdTrapezoidalArea(double, double, double);
 void PropertiesAuxiliary(BODY *, CONTROL *, SYSTEM *, UPDATE *);
 void fdGetUpdateInfo(BODY *, CONTROL *, SYSTEM *, UPDATE *,
                      fnUpdateVariable ***);

src/options.c

@@ -1078,7 +1078,6 @@ void ReadUnitTemp(CONTROL *control, FILES *files, OPTIONS *options, int iFile) {
 
 void ReadSystemName(CONTROL *control, FILES *files, OPTIONS *options,
                     SYSTEM *system, int iFile) {
-  /* System Name */
   int lTmp = -1;
   char cTmp[OPTLEN];
 
@@ -1092,6 +1091,28 @@ void ReadSystemName(CONTROL *control, FILES *files, OPTIONS *options,
   }
 }
 
+void ReadSystemAge(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
+                   SYSTEM *system, int iFile) {
+  int lTmp = -1;
+  double dTmp;
+
+  AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
+                  control->Io.iVerbose);
+  if (lTmp >= 0) {
+    CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp,
+                     control->Io.iVerbose);
+    if (dTmp < 0) {
+      system->dAge = dTmp * dNegativeDouble(*options, files->Infile[iFile].cIn,
+                                            control->Io.iVerbose);
+    } else {
+      system->dAge = dTmp * fdUnitsTime(control->Units[iFile].iTime);
+    }
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
+  } else if (iFile > 0) {
+    AssignDefaultDouble(options, &system->dAge, files->iNumInputs);
+  }
+}
+
 void ReadBodyFileNames(BODY **body, CONTROL *control, FILES *files,
                        OPTIONS *options, char *cFile, char ***saBodyFiles,
                        int *iNumLines, int *iaLines) {
@@ -1442,8 +1463,6 @@ void ReadOutputTime(BODY *body, CONTROL *control, FILES *files,
   }
 }
 
-/* Backward integration stop time */
-
 void ReadStopTime(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
                   SYSTEM *system, int iFile) {
   /* This parameter can exist in any file, but only once */
@@ -1470,6 +1489,32 @@ void ReadStopTime(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
   }
 }
 
+void ReadStopAge(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
+                 SYSTEM *system, int iFile) {
+  /* This parameter can exist in any file, but only once */
+  int lTmp = -1;
+  double dTmp;
+
+  AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
+                  control->Io.iVerbose);
+  if (lTmp >= 0) {
+    /* Option was found */
+    CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp,
+                     control->Io.iVerbose);
+    if (dTmp < 0) {
+      if (control->Io.iVerbose >= VERBERR) {
+        fprintf(stderr, "ERROR: %s must be greater than 0.\n", options->cName);
+      }
+      LineExit(files->Infile[iFile].cIn, lTmp);
+    }
+    /* Convert stop time to cgs */
+    control->Evolve.dStopAge = dTmp * fdUnitsTime(control->Units[iFile].iTime);
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
+  } else {
+    AssignDefaultDouble(options, &control->Evolve.dStopAge, files->iNumInputs);
+  }
+}
+
 /* Integration timestep */
 
 void ReadTimeStep(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
@@ -2617,7 +2662,7 @@ void fvAllocateOutputArrays(char ****saMatch, char ***saOutput, int **baNeg,
 
   for (iIndex = 0; iIndex < iNumArgs; iIndex++) {
     (*saMatch)[iIndex] =
-          malloc(MAXARRAY * sizeof(char*)); // Could be this many matches
+          malloc(MAXARRAY * sizeof(char *)); // Could be this many matches
     for (iMatch = 0; iMatch < MAXARRAY; iMatch++) {
       (*saMatch)[iIndex][iMatch] = NULL;
     }
@@ -3712,7 +3757,7 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
    */
 
   fvFormattedString(&options[OPT_AGE].cName, "dAge");
-  fvFormattedString(&options[OPT_AGE].cDescr, "System Age");
+  fvFormattedString(&options[OPT_AGE].cDescr, "Body's Age");
   fvFormattedString(&options[OPT_AGE].cDefault, "0");
   fvFormattedString(&options[OPT_AGE].cNeg, "Gyr");
   fvFormattedString(&options[OPT_AGE].cDimension, "time");
@@ -3724,6 +3769,20 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_AGE].iFileType  = 2;
   fnRead[OPT_AGE]             = &ReadAge;
 
+  int iOpt = OPT_SYSTEMAGE;
+  fvFormattedString(&options[iOpt].cName, "dSystemAge");
+  fvFormattedString(&options[iOpt].cDescr, "System Age");
+  fvFormattedString(&options[iOpt].cDefault, "0");
+  fvFormattedString(&options[iOpt].cNeg, "1 Gyr");
+  fvFormattedString(&options[iOpt].cDimension, "time");
+  options[iOpt].dDefault   = 0;
+  options[iOpt].iType      = 2;
+  options[iOpt].iModuleBit = 0;
+  options[iOpt].bNeg       = 1;
+  options[iOpt].dNeg       = 1e9 * YEARSEC;
+  options[iOpt].iFileType  = 2;
+  fnRead[iOpt]             = &ReadSystemAge;
+
   fvFormattedString(&options[OPT_ALBEDOGLOBAL].cName, "dAlbedoGlobal");
   fvFormattedString(&options[OPT_ALBEDOGLOBAL].cDescr,
                     "Globally averaged albedo");
@@ -3805,6 +3864,20 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_STOPTIME].iFileType  = 2;
   fnRead[OPT_STOPTIME]             = &ReadStopTime;
 
+  iOpt = OPT_STOPAGE;
+  fvFormattedString(&options[iOpt].cName, "dStopAge");
+  fvFormattedString(&options[iOpt].cDescr, "Age to stop integration");
+  fvFormattedString(&options[iOpt].cDefault, "10 Gigayears");
+  fvFormattedString(&options[iOpt].cNeg, "Years");
+  fvFormattedString(&options[iOpt].cDimension, "time");
+  options[iOpt].dDefault   = 1e10 * YEARSEC;
+  options[iOpt].iType      = 2;
+  options[iOpt].iModuleBit = 0;
+  options[iOpt].bNeg       = 1;
+  options[iOpt].dNeg       = YEARSEC;
+  options[iOpt].iFileType  = 2;
+  fnRead[iOpt]             = &ReadStopAge;
+
   fvFormattedString(&options[OPT_TIMESTEP].cName, "dTimeStep");
   fvFormattedString(&options[OPT_TIMESTEP].cDescr, "Integration Timestep");
   fvFormattedString(&options[OPT_TIMESTEP].cDefault, "1 year");

src/options.h

@@ -30,13 +30,15 @@
 // Regular Options
 
 #define OPT_AGE 100
+#define OPT_SYSTEMAGE 102
 #define OPT_ALBEDOGLOBAL 105
 
 #define OPT_BACK 110
 #define OPT_OUTFILE 120
 #define OPT_ETA 130
 #define OPT_OUTPUTTIME 140
 #define OPT_STOPTIME 150
+#define OPT_STOPAGE 155
 #define OPT_TIMESTEP 160
 #define OPT_VARDT 170
 #define OPT_BODYNAME 180
@@ -128,8 +130,7 @@ void ReadOptions(BODY **, CONTROL *, FILES *, MODULE *, OPTIONS *, OUTPUT *,
                  SYSTEM *, UPDATE **, fnReadOption *, char[]);
 
 double dNegativeDouble(OPTIONS, char[], int);
-void AddOptionStringArray(char[], char[], char***, int *, int *,
-                          int *, int);
+void AddOptionStringArray(char[], char[], char ***, int *, int *, int *, int);
 void AddOptionDoubleArray(char[], char[], double *, int *, int *, int *, int);
 void NotPrimaryInput(int, char[], char[], int, int);
 void AddOptionDouble(char[], char[], double *, int *, int);