diff --git a/src/distrot.c b/src/distrot.c index c263d979c..a42e524dd 100644 --- 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; } /**