Merge pull request #275 from deitrr/main

Precession and memory fixes in DistRot for ReadOrbitData use
VirtualPlanetaryLaboratory · Mar 2, 2024 · 76f61f8 · 76f61f8
2 parents cc7a60a + 69e6ab5
commit 76f61f8
Showing 1 changed file with 35 additions and 15 deletions.
diff --git a/src/distrot.c b/src/distrot.c
@@ -420,9 +420,9 @@ void VerifyOrbitData(BODY *body, CONTROL *control, OPTIONS *options,
       }
       // Check file has exactly 7 columns
       if (fgets(cLine, LINE, fileorb) == NULL) {
-        fprintf(stderr,"ERROR: Unable to read line from orbit data file.");
+        fprintf(stderr, "ERROR: Unable to read line from orbit data file.");
         exit(EXIT_INPUT);
-      } 
+      }
       GetWords(cLine, cFoo, &iNumColsFound, &bFoo);
       if (iNumCols != iNumColsFound) {
         if (control->Io.iVerbose >= VERBERR) {
@@ -459,10 +459,10 @@ void VerifyOrbitData(BODY *body, CONTROL *control, OPTIONS *options,
 
       iLine = 0;
       while (feof(fileorb) == 0) {
-        if (fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf\n", &dttmp, &datmp, &detmp,
-               &ditmp, &daptmp, &dlatmp, &dmatmp) != 7) {
-                  fprintf(stderr,"ERROR: Incorrect number of columns in orbit file.");
-                  exit(EXIT_INPUT);
+        if (fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf\n", &dttmp, &datmp,
+                   &detmp, &ditmp, &daptmp, &dlatmp, &dmatmp) != 7) {
+          fprintf(stderr, "ERROR: Incorrect number of columns in orbit file.");
+          exit(EXIT_INPUT);
         }
         body[iBody].daTimeSeries[iLine] =
               dttmp * fdUnitsTime(control->Units[iBody + 1].iTime);
@@ -907,12 +907,20 @@ void FinalizeUpdateZoblDistRot(BODY *body, UPDATE *update, int *iEqn, int iVar,
                                int iBody, int iFoo) {
   int iPert;
 
-  update[iBody].padDZoblDtDistRot =
-        malloc((body[iBody].iGravPerts) * sizeof(double *));
-  update[iBody].iaZoblDistRot = malloc((body[iBody].iGravPerts) * sizeof(int));
-  for (iPert = 0; iPert < body[iBody].iGravPerts; iPert++) {
+  if (body[iBody].bReadOrbitData) {
+    update[iBody].padDZoblDtDistRot     = malloc(1 * sizeof(double *));
+    update[iBody].iaZoblDistRot         = malloc(1 * sizeof(int));
     update[iBody].iaModule[iVar][*iEqn] = DISTROT;
-    update[iBody].iaZoblDistRot[iPert]  = (*iEqn)++;
+    update[iBody].iaZoblDistRot[0]      = (*iEqn)++;
+  } else {
+    update[iBody].padDZoblDtDistRot =
+          malloc((body[iBody].iGravPerts) * sizeof(double *));
+    update[iBody].iaZoblDistRot =
+          malloc((body[iBody].iGravPerts) * sizeof(int));
+    for (iPert = 0; iPert < body[iBody].iGravPerts; iPert++) {
+      update[iBody].iaModule[iVar][*iEqn] = DISTROT;
+      update[iBody].iaZoblDistRot[iPert]  = (*iEqn)++;
+    }
   }
 }
 
@@ -2110,13 +2118,19 @@ external model
 @return Derivative dx/dt
 */
 double fndDistRotExtDxDt(BODY *body, SYSTEM *system, int *iaBody) {
-  double y;
+  double y, dprec;
   y = fabs(1.0 - (body[iaBody[0]].dXobl * body[iaBody[0]].dXobl) -
            (body[iaBody[0]].dYobl * body[iaBody[0]].dYobl));
 
+  if (body[iaBody[0]].bForcePrecRate == 0) {
+    dprec = fndCentralTorqueR(body, iaBody[0]);
+  } else {
+    dprec = body[iaBody[0]].dPrecRate;
+  }
+
   return fndObliquityAExt(body, system, iaBody) * sqrt(y) +
          body[iaBody[0]].dYobl * 2. * fndObliquityCExt(body, system, iaBody) -
-         body[iaBody[0]].dYobl * fndCentralTorqueR(body, iaBody[0]);
+         body[iaBody[0]].dYobl * dprec;
 }
 
 /**
@@ -2128,13 +2142,19 @@ external model
 @return Derivative dy/dt
 */
 double fndDistRotExtDyDt(BODY *body, SYSTEM *system, int *iaBody) {
-  double y;
+  double y, dprec;
   y = fabs(1.0 - (body[iaBody[0]].dXobl * body[iaBody[0]].dXobl) -
            (body[iaBody[0]].dYobl * body[iaBody[0]].dYobl));
 
+  if (body[iaBody[0]].bForcePrecRate == 0) {
+    dprec = fndCentralTorqueR(body, iaBody[0]);
+  } else {
+    dprec = body[iaBody[0]].dPrecRate;
+  }
+
   return -fndObliquityBExt(body, system, iaBody) * sqrt(y) -
          body[iaBody[0]].dXobl * 2. * fndObliquityCExt(body, system, iaBody) +
-         body[iaBody[0]].dXobl * fndCentralTorqueR(body, iaBody[0]);
+         body[iaBody[0]].dXobl * dprec;
 }
 
 /**