sStellarModel SineWave is working with optimization, but the luminosi…

…ty becomes negative when running in debug mode.
VirtualPlanetaryLaboratory · Feb 13, 2024 · 4b89543 · 4b89543
1 parent f45daf1
commit 4b89543
Show file tree

Hide file tree

Showing 7 changed files with 27 additions and 11 deletions.
diff --git a/examples/LuminosityCycle/none.in b/examples/LuminosityCycle/none.in
@@ -3,5 +3,6 @@ sName	                  none
 saModules	              stellar
 sStellarModel             none
 dMass                     1
+dRadius                   0.00454
 dLuminosity               -1
 saOutputOrder Time -Luminosity -Radius Temperature -RotPer -LXUVTot RadGyra
diff --git a/examples/LuminosityCycle/sinewave.in b/examples/LuminosityCycle/sinewave.in
@@ -3,9 +3,10 @@ sName	                  sinewave
 saModules	              stellar
 sStellarModel             sinewave
 dMass                     1
+dRadius                   0.00454
 dLuminosity               -1
 dLuminosityAmplitude      -0.1
-dLuminosityPeriod         -10
-dLuminosityPhase          30
+dLuminosityPeriod         -100
+dLuminosityPhase          0
 
 saOutputOrder Time -Luminosity -Radius Temperature -RotPer -LXUVTot RadGyra
diff --git a/examples/LuminosityCycle/vpl.in b/examples/LuminosityCycle/vpl.in
@@ -13,5 +13,5 @@ dMinValue                 1e-10
 bDoForward                1
 bVarDt                    1
 dEta                      0.01
-dStopTime                 5e9
-dOutputTime               1e6
+dStopTime                 1000
+dOutputTime               1
diff --git a/src/body.c b/src/body.c
@@ -461,6 +461,7 @@ void BodyCopy(BODY *dest, BODY *src, EVOLVE *evolve) {
     dest[iBody].dLostEng      = src[iBody].dLostEng;
     dest[iBody].dAlbedoGlobal = src[iBody].dAlbedoGlobal;
     dest[iBody].bCalcDynEllip = src[iBody].bCalcDynEllip;
+    dest[iBody].dRadGyra = src[iBody].dRadGyra;
 
     dest[iBody].bBinary   = src[iBody].bBinary;
     dest[iBody].bDistOrb  = src[iBody].bDistOrb;
@@ -1488,3 +1489,8 @@ void fdHabitableZoneKopparapu2013(BODY *body, int iNumBodies,
     daHZLimit[iLimit] = pow(dLuminosity / dSeff[iLimit], 0.5) * AUM;
   }
 }
+
+double fdEffectiveTemperature(BODY *body,int iBody) {
+  double dTeff = pow((body[iBody].dLuminosity/(4*PI*SIGMA*body[iBody].dRadius*body[iBody].dRadius)),0.25);
+  return dTeff;
+}
diff --git a/src/body.h b/src/body.h
@@ -108,6 +108,8 @@ double CalcDynEllipEq(BODY *, int);
 
 void fdHabitableZoneKopparapu2013(BODY *, int, double *);
 
+double fdEffectiveTemperature(BODY*,int);
+
 // RB: Move
 
 // Proxima Centauri properties

diff --git a/src/stellar.c b/src/stellar.c
@@ -765,7 +765,7 @@ void VerifyStellarSineWave(BODY *body, CONTROL *control, OPTIONS *options,
              options[OPT_STELLARMODEL].iLine[iBody + 1]);
   }
 
-  body[iBody].dLuminosity = fdLuminosityFunctionSineWave(body, iBody);
+  body[iBody].dLuminosityInitial = body[iBody].dLuminosity;
 }
 
 void VerifyStellarNone(BODY *body, CONTROL *control, OPTIONS *options,
@@ -1496,7 +1496,8 @@ double fdRadius(BODY *body, SYSTEM *system, int *iaBody) {
     }
   }
   if (body[iaBody[0]].iStellarModel == STELLAR_MODEL_NONE ||
-      body[iaBody[0]].iStellarModel == STELLAR_MODEL_CONST) {
+      body[iaBody[0]].iStellarModel == STELLAR_MODEL_CONST ||
+      body[iaBody[0]].iStellarModel == STELLAR_MODEL_SINEWAVE) {
     return body[iaBody[0]].dRadius;
   } else {
     return 0;
@@ -1522,6 +1523,10 @@ double fdTemperature(BODY *body, SYSTEM *system, int *iaBody) {
       body[iaBody[0]].iStellarModel = STELLAR_MODEL_CONST;
     }
   }
+  if (body[iaBody[0]].iStellarModel == STELLAR_MODEL_SINEWAVE) {
+    foo = fdEffectiveTemperature(body,iaBody[0]);
+    return foo;
+  }
   if (body[iaBody[0]].iStellarModel == STELLAR_MODEL_NONE ||
       body[iaBody[0]].iStellarModel == STELLAR_MODEL_CONST) {
     return body[iaBody[0]].dTemperature;
@@ -1554,7 +1559,8 @@ double fdRadGyra(BODY *body, SYSTEM *system, int *iaBody) {
     }
   }
   if (body[iaBody[0]].iStellarModel == STELLAR_MODEL_NONE ||
-      body[iaBody[0]].iStellarModel == STELLAR_MODEL_CONST) {
+      body[iaBody[0]].iStellarModel == STELLAR_MODEL_CONST||
+      body[iaBody[0]].iStellarModel == STELLAR_MODEL_SINEWAVE) {
     return body[iaBody[0]].dRadGyra;
   } else {
     return 0;
@@ -2012,9 +2018,8 @@ double fdSurfEnFluxStellar(BODY *body, SYSTEM *system, UPDATE *update,
 
 double fdLuminosityFunctionSineWave(BODY *body, int iBody) {
   double dLuminosity =
-        body[iBody].dLuminosity +
-        body[iBody].dLuminosityAmplitude *
-              sin(body[iBody].dAge * body[iBody].dLuminosityFrequency +
-                  body[iBody].dLuminosityPhase);
+      body[iBody].dLuminosityInitial + body[iBody].dLuminosityAmplitude *
+      sin(body[iBody].dAge * body[iBody].dLuminosityFrequency +
+      body[iBody].dLuminosityPhase);
   return dLuminosity;
 }
diff --git a/src/vplanet.h b/src/vplanet.h
@@ -614,6 +614,7 @@ struct BODY {
                          Ro>ROSSBYCRIT */
   int bEvolveRG; /**< Whether or not to evolve radius of gyration? Defaults to 0
                   */
+  double dLuminosityInitial;           
   double dLuminosityAmplitude;
   double dLuminosityFrequency;
   double dLuminosityPhase;