Merge pull request #26 from VirtualPlanetaryLaboratory/dev

Dev
VirtualPlanetaryLaboratory · Dec 29, 2018 · 7a72c80 · 7a72c80
2 parents 009138e + 7d41da1
commit 7a72c80
Show file tree

Hide file tree

Showing 7 changed files with 44 additions and 3 deletions.
diff --git a/src/module.c b/src/module.c
@@ -1267,11 +1267,13 @@ void  VerifyModuleCompatability(BODY *body,CONTROL *control,FILES *files,MODULE
         fprintf(stderr,"ERROR: Modules Stellar and RadHeat cannot be applied to the same body.\n");
       LineExit(files->Infile[iBody+1].cIn,options[OPT_MODULES].iLine[iBody+1]);
     }
+    /* This combo seems OK, and is necessary to calculate instellation. I'll leave it
+        commented for now.
     if (body[iBody].bSpiNBody) {
       if (control->Io.iVerbose >= VERBERR)
         fprintf(stderr,"ERROR: Modules Stellar and SpiNBody cannot be applied to the same body.\n");
       LineExit(files->Infile[iBody+1].cIn,options[OPT_MODULES].iLine[iBody+1]);
-    }
+    } */
     if (body[iBody].bThermint) {
       if (control->Io.iVerbose >= VERBERR)
         fprintf(stderr,"ERROR: Modules Stellar and ThermInt cannot be applied to the same body.\n");

diff --git a/src/options.h b/src/options.h
@@ -93,7 +93,6 @@
 #define OPT_PRILUM	        615
 #define OPT_SPECMOMINERTIA       616
 
-
 #define OPT_VISCUMAN            630   //Viscosity UMTBL
 
 #define OPT_RADIUS              650

diff --git a/src/output.c b/src/output.c
@@ -124,6 +124,24 @@ void WriteHZLimitDryRunaway(BODY *body,CONTROL *control,OUTPUT *output,SYSTEM *s
    }
  }
 
+ void WriteInstellation(BODY *body,CONTROL *control,OUTPUT *output,SYSTEM *system,UNITS *units,UPDATE *update,int iBody,double *dTmp,char cUnit[]) {
+
+   // Should have special case if bBinary=1
+   if (body[iBody].bSpiNBody)
+    //Bary2OrbElems(body,control->Evolve.iNumBodies);
+
+   *dTmp = fdInstellation(body,iBody);
+
+   if (output->bDoNeg[iBody]) {
+     *dTmp *= output->dNeg;
+     strcpy(cUnit,output->cNeg);
+   } else {
+     *dTmp /= fdUnitsAngle(units->iAngle);
+     fsUnitsAngle(units->iAngle,cUnit);
+   }
+ }
+
+
 /*
  * K
  */
@@ -897,6 +915,15 @@ void InitializeOutputGeneral(OUTPUT *output,fnWriteOutput fnWrite[]) {
   output[OUT_INC].iModuleBit = DISTORB + GALHABIT + SPINBODY + BINARY;
   fnWrite[OUT_INC] = &WriteBodyInc;
 
+  sprintf(output[OUT_INSTELLATION].cName,"Instellation");
+  sprintf(output[OUT_INSTELLATION].cDescr,"Orbit-averaged INcident STELLar radiATION");
+  sprintf(output[OUT_INSTELLATION].cNeg,"W/m^2");
+  output[OUT_INSTELLATION].bNeg = 1;
+  output[OUT_INSTELLATION].dNeg = 1;
+  output[OUT_INSTELLATION].iNum = 1;
+  output[OUT_INSTELLATION].iModuleBit = 1;
+  fnWrite[OUT_INSTELLATION] = &WriteInstellation;
+
   /*
    * K
    */

diff --git a/src/output.h b/src/output.h
@@ -62,6 +62,8 @@
 #define OUT_HECC                600
 #define OUT_HZLIMITDRYRUNAWAY   602
 
+#define OUT_INSTELLATION        604
+
 #define OUT_KECC                605
 #define OUT_ORBECC		          610
 #define OUT_ORBEN		            620

diff --git a/src/system.c b/src/system.c
@@ -295,6 +295,15 @@ void CalcPQ(BODY *body, int iBody) {
   body[iBody].dQinc = body[iBody].dSinc*cos(body[iBody].dLongA);
 }
 
+/* Instellation */
+double fdInstellation(BODY *body,int iBody) {
+
+  if (iBody > 0)
+    return body[0].dLuminosity/(4*PI*body[iBody].dSemi*body[iBody].dSemi*sqrt(1-body[iBody].dEcc*body[iBody].dEcc));
+  else // Central body can't have instellation (for now) XXX
+  return -1;
+}
+
 /*
  *
  * Multi-body matrix equations

diff --git a/src/system.h b/src/system.h
@@ -29,6 +29,8 @@ int bPrimary(BODY*,int);
 void CalcHK(BODY*,int);
 void CalcPQ(BODY*,int);
 
+double fdInstellation(BODY*,int);
+
 /* Multi-body derivative equations */
 double fdSemiDtEqSt(BODY*,SYSTEM*,int*);
 double fdJStarDt(BODY*,int);

diff --git a/src/vplanet.h b/src/vplanet.h
@@ -70,11 +70,11 @@
 #define DAYSEC        86400       // Seconds per day
 #define REARTH        6.3781e6    // Equatorial; Prsa et al. 2016
 #define RJUP          7.1492e7    // Equatorial; Prsa et al. 2016
-#define MJUP          1.898130e27 // Prsa et al. 2016
 #define RNEP          2.4764e7    // Neptune's Radius (ref?)
 #define MNEP          1.0244e26   // Neptune's Mass (ref?)
 #define RHOEARTH      5515        // Earth's Density
 #define eEARTH        0.016710219 // Earth's Eccentricity
+#define MJUP          1.898130e27 // Prsa et al. 2016
 #define YEARDAY       365.25      // Days per year (more precise??)
 #define MSAT          5.6851e26   // Saturns' Mass (ref?)
 #define DEGRAD        0.017453292519444445 // Degrees per radian