EQTIDE
eqtide.h
1 /**********************************************************************
2  * *
3  * EQTIDE is a general integrator of "equilibrium tide" models *
4  * as originally conceived by George Darwin. These equation were *
5  * first assembled in Heller, Leconte & Barnes (2011), and the *
6  * version first is released simultaneously with Barnes, R. (2016). *
7  * *
8  * EQTIDE is optimized for a rocky planet orbiting a main sequence *
9  * star, but is easily applied to any generic two-body system. *
10  * Example runs are provided at https://github.com/RoryBarnes/EqTide. *
11  * Please submit pull requests if you find an error! *
12  * *
13  * EQTIDE was written by Rory Barnes. *
14  * *
15  * eqtide.h *
16  * *
17  * Header file for EQTIDE. All structs and global variables are *
18  * declared here and in options.h and output.h. *
19  * *
20  **********************************************************************/
21 
22 #define MEARTH 5.972186e27 // Prsa et al. 2016
23 #define MSUN 1.988416e33 // Prsa et al. 2016
24 #define AUCM 1.49598e13
25 #define RSUN 6.957e10 // Prsa et al. 2016
26 #define YEARSEC 3.15576e7
27 #define DAYSEC 86400
28 #define REARTH 6.3781e8 // Equatorial; Prsa et al. 2016
29 #define RJUP 7.1492e9 // Equatorial; Prsa et al. 2016
30 #define MJUP 1.898130e30 // Prsa et al. 2016
31 
32 #define BIGG 6.67428e-8 // From Luzum et al., 2011; value recommended by IAU NSFA in Prsa et al. 2016
33 #define PI 3.1415926535
34 
35 #define RNEP 2.4764e9
36 #define MNEP 1.0244e29
37 #define RHOEARTH 5.52
38 #define YEARDAY 365.25
39 #define MSAT 5.6851e29
40 #define DEGRAD 0.017453292519444445
41 
42 #define EXIT_EXE 1
43 #define EXIT_PARAM 2
44 #define EXIT_UNITS 3
45 #define EXIT_WRITE 4
46 #define EXIT_INT 5
47 #define EXIT_OUTPUT 6
48 
49 #define VERBERR 1
50 #define VERBPROG 2
51 #define VERBINPUT 3
52 #define VERBUNITS 4
53 #define VERBALL 5
54 
55 #define OPTLEN 24
56 #define LINELEN 256
57 #define NUMOUT 1000 /* Number of output parameters */
58 #define NUMOPT 1000
59 
60 /* EPS is the floor for which two numbers are deemed to be unequal. If the
61 difference between two numbers is < EPS, then they are assumed to be the same.
62 This is particularly important for tidally locked body for which round-off
63 errors can cause the difference of the mean motion and rotational frequency
64 to be non-zero. */
65 #define EPS 1e-10
66 
67 #define TINY (1./HUGE)
68 
69 #define CPL 1
70 #define CTL 10
71 
72 #define EULER 0
73 #define RK4 1
74 
78 typedef struct {
79  double dMass;
80  double dRadius;
81  double dTau;
82  double dQ;
83  double dK2;
84  double dObliquity;
85  double dSpinRate;
86  double dRG;
87  double *iEpsilon;
89  int bForceEqSpin;
90  double dMaxLockDiff;
92  int iMassRad;
94  double dDomegaDt;
95  double dDobliquityDt;
96 } PRIMARY;
97 
101 typedef struct {
102  double dMass;
103  double dRadius;
104  double dTau;
105  double dQ;
106  double dK2;
107  double dObliquity;
108  double dSpinRate;
109  double dRG;
110  double dSyncEcc;
112  int bForceEqSpin;
113  double dMaxLockDiff;
115  double dSemi;
116  double dEcc;
117  double dMeanMotion;
118  double dAge;
120  int iMassRad;
122  double dDomegaDt;
123  double dDobliquityDt;
124  double dDaDt;
125  double dDeDt;
126 } SECONDARY;
127 
128 typedef struct {
129  int iVerbose;
131  int iDigits;
132  int iSciNot;
134  int exit_param;
135  int exit_io;
136  int exit_input;
137  int exit_exe;
138  int exit_units;
139  int exit_output;
141 } IO;
142 
146 typedef void (*fvDerivs)(PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double,int);
147 
150 typedef double (*fdEqSpin)(double,double,double,int);
151 
152 typedef struct {
153  int bHalt;
154  int bDblSync;
155  int bMerge;
156  double dMinSemi;
157  double dMinPriObl;
158  double dMinSecObl;
159  double dMaxEcc;
160  double dMinEcc;
161  int bPosDeDt;
162  int bPriLock;
163  int bSecLock;
164  int bSecSync;
165 } HALT;
166 
167 typedef struct {
168  HALT halt;
169 
170  int iUnitMass;
171  int iUnitLength;
172  int iUnitAngle;
173  int iUnitTime;
175  char cSystemName[16];
176  int iTideModel;
177  int bDiscreteRot;
179  int bLog;
181  double l0;
182  double dl;
184  char *cOutputOrder[NUMOUT];
185  int iNumCols;
187  double dMinValue;
189  /* Parameters to perform forward integration */
190  int bDoForward;
191  double dForwTimeStep;
192  double dForwStopTime;
193  double dForwOutputTime;
195  /* Parameters to perform reverse integration */
196  int bDoBackward;
197  double dBackTimeStep;
198  double dBackStopTime;
199  double dBackOutputTime;
201  int bVarDt;
202  double dTimestepCoeff;
204  int iIntegration;
205  fvDerivs fDerivs;
206  fdEqSpin fEqSpin;
207 } PARAM;
208 
209 
211 typedef double (*fdStep)(PARAM*,PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double*,double,int);
212 
213 typedef struct {
214  char *cExe;
215  char *cIn;
216  char cFor[32];
217  char cBack[32];
218  char cLog[32];
219 } FILES;
220 
221 typedef struct {
222  char cParam[NUMOPT][OPTLEN];
223  char cDescr[NUMOPT][LINELEN];
224  int iType[NUMOPT];
225  char cDefault[NUMOPT][OPTLEN];
226 } OPTIONS;
227 
228 typedef struct {
229  char *cParam[NUMOUT];
230  char cDescr[NUMOUT][LINELEN];
231  int iNeg[NUMOUT];
232  char cNeg[NUMOUT][OPTLEN];
233  double dConvert[NUMOUT];
234 } OUTPUT;
235 
236 // String/print utilities
237 char *lower(char[]);
238 void fprintd(FILE*,double,int,int);
239 char *InitializeString();
240 
241 /* Mass-Radius relations */
242 double dBaylessOrosz06_MassRad(double);
243 double dGordaSvech99_MassRad(double);
244 double dReidHawley_MassRad(double);
245 double dSotin07_MassRad(double);
246 double dBaylessOrosz06_RadMass(double);
247 double dBaraffe15_MassRad(double);
248 double dGordaSvech99_RadMass(double);
249 double dReidHawley_RadMass(double);
250 double dSotin07_RadMass(double);
251 
252 double dDpDt(double,double,double);
253 void AssignChi(PRIMARY*,SECONDARY*,double*);
254 double a2p(double,double);
255 double p2a(double,double);
256 double dFreqToPer(double);
257 double dPerToFreq(double);
258 double dSemiToMeanMotion(double,double);
259 double dRotVel(double,double);
260 double DOrbPerDt(double,double,double);
261 
262 /* Units */
263 double dLengthUnit(int,int);
264 double dTimeUnit(int,int);
265 double dMassUnit(int,int);
266 double dAngleUnit(int,int);
267 
268 /* CPL2 Functions */
269 double EqSpinRate_CPL2(double,double,double,int);
270 void AssignEpsilon(double,double,int*);
271 void AssignZprime(PRIMARY*,SECONDARY*,double *);
272 double dDaDtAnn_CPL2(PRIMARY*,SECONDARY*,int**,double *);
273 double dDaDt_CPL2(PRIMARY*,SECONDARY*,int **,double*);
274 double dDeDt_CPL2(PRIMARY*,SECONDARY*,int **,double*);
275 double dDomegaDt_CPL2(double,double,double,double,double,double,int*,double);
276 double dDoblDt_CPL2(double,double,double,double,double,int*,double,double,double);
277 double dTideHeat_CPL2(int*,double,double,double,double,double);
278 void TideRewind_CPL2(PARAM*,PRIMARY*,SECONDARY*,OUTPUT*,FILES*,IO*,fdStep);
279 void TideForward_CPL2(PARAM*,PRIMARY*,SECONDARY*,OUTPUT*,FILES*,IO*,fdStep);
280 double dDaDt1_CPL2(double,double,double,double,double,int*,double);
281 double dDeDt1_CPL2(double,double,double,double,int*,double);
282 double dTideHeatEq_CPL2(double,double,double,double,int);
283 double EqSpinRate_CPL2Discrete(double,double);
284 double EqSpinRate_CPL2Cont(double,double);
285 double dGammaRot(double,double,int*);
286 double dGammaOrb(double,double,int*);
287 void DerivsCPL(PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double,int);
288 
289 /* CTL8 Functions */
290 double EqSpinRate_CTL8(double,double,double,int);
291 double AssignF1(double);
292 double AssignF2(double);
293 double AssignF3(double);
294 double AssignF4(double);
295 double AssignBeta(double);
296 double AssignF5(double);
297 void AssignZ(PRIMARY*,SECONDARY*,double *);
298 double dDaDt_CTL8(PRIMARY*,SECONDARY*,double *,double,double *);
299 double dDeDt_CTL8(PRIMARY*,SECONDARY*,double *,double,double *);
300 double dDomegaDt_CTL8(double,double,double,double,double,double,double,double,double,double *);
301 double dDoblDt_CTL8(double,double,double,double,double,double,double,double,double,double*,double);
302 double dAnnualHeat_CTL8(double,double,double);
303 double dDailyHeat_CTL8(double,double*,double,double,double,double);
304 void TideForward_CTL8 (PARAM*,PRIMARY*,SECONDARY*,OUTPUT*,FILES*,IO*,fdStep);
305 void TideRewind_CTL8 (PARAM*,PRIMARY*,SECONDARY*,OUTPUT*,FILES*,IO*,fdStep);
306 double dDaDt1_CTL8(double,double,double,double,double,double,double,double,double *,double);
307 double dDeDt1_CTL8(double,double,double,double,double,double,double,double,double*,double);
308 double dTideHeat_CTL8(double,double*,double,double,double,double);
309 double dTideHeatEq_CTL8(double,double*,double,double,double);
310 void DerivsCTL(PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double,int);
311 
312 /* Integration Methods */
313 double EulerStep(PARAM*,PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double*,double,int);
314 double RK4Step(PARAM*,PRIMARY*,SECONDARY*,IO*,double*,double*,double*,double,int**,double*,double,int);
315 
316 /* Conserved Quantities */
317 double dOrbAngMom(PRIMARY*,SECONDARY*);
318 double dSpinAngMom(double,double,double,double);
319 double dRotEn(double,double,double,double);
320 double dOrbEn(double,double,double);
321 
322 void LineExit(char[],int, int, int);
323 int iSign(double);
324 double dDPerDt(double,double);
325 
326 // Major Functions
327 void InitializeIO(IO*);
328 void InitializeOutput(OUTPUT*);
329 void InitializeOptions(OPTIONS*);
330 void WriteHelp(OUTPUT,OPTIONS,char[]);
331 void ReadOptions(char[],OPTIONS,PARAM*,PRIMARY*,SECONDARY*,FILES*,OUTPUT*,IO*,fdStep*);
332 void WriteLog(PARAM,PRIMARY,SECONDARY,OUTPUT,FILES,IO,int);
333 void Output(PARAM*,PRIMARY*,SECONDARY*,OUTPUT*,IO*,double,double,FILE*);
PRIMARY::dDobliquityDt
double dDobliquityDt
Time derivative of obliquity.
Definition: eqtide.h:95
IO::exit_exe
int exit_exe
Exit code for execution error.
Definition: eqtide.h:137
SECONDARY::dDaDt
double dDaDt
Time derivative of semi-major axis.
Definition: eqtide.h:124
OUTPUT
Definition: eqtide.h:228
IO
Definition: eqtide.h:128
SECONDARY::dRG
double dRG
Radius of Gyration.
Definition: eqtide.h:109
SECONDARY::dMass
double dMass
Planetary Mass.
Definition: eqtide.h:102
SECONDARY::dSyncEcc
double dSyncEcc
Synchronous rotation below this e.
Definition: eqtide.h:110
PRIMARY::dK2
double dK2
Love number of primary.
Definition: eqtide.h:83
IO::iVerbose
int iVerbose
Verbosity: 0=no STDOUT, 1=all.
Definition: eqtide.h:129
SECONDARY::dAge
double dAge
System Age.
Definition: eqtide.h:118
SECONDARY::dDobliquityDt
double dDobliquityDt
Time derivative of obliquity.
Definition: eqtide.h:123
HALT::bSecSync
int bSecSync
Halt if secondary rotation becomes synchronous?
Definition: eqtide.h:164
SECONDARY::iMassRad
int iMassRad
Mass-Radius relationship identifier.
Definition: eqtide.h:120
PARAM::bDiscreteRot
int bDiscreteRot
Use Discrete Rotation model (CPL)
Definition: eqtide.h:177
SECONDARY::dQ
double dQ
Secondary Tidal Q.
Definition: eqtide.h:105
IO::exit_param
int exit_param
Exit code for failure to specify necessary paramter.
Definition: eqtide.h:134
PARAM
Definition: eqtide.h:167
SECONDARY::bForceEqSpin
int bForceEqSpin
Tidally lock second?
Definition: eqtide.h:112
PRIMARY::dDomegaDt
double dDomegaDt
Time derivative of spin rate.
Definition: eqtide.h:94
HALT::bDblSync
int bDblSync
Halt for double synchronous rotation?
Definition: eqtide.h:154
PARAM::dBackOutputTime
double dBackOutputTime
Output interval in Backward integration (years)
Definition: eqtide.h:199
PARAM::dBackStopTime
double dBackStopTime
Backward integration time (years)
Definition: eqtide.h:198
IO::exit_units
int exit_units
Exit code for unit conversion error.
Definition: eqtide.h:138
IO::iSciNot
int iSciNot
Decade to change from flt.
Definition: eqtide.h:132
SECONDARY::dRadius
double dRadius
Secondary Radius.
Definition: eqtide.h:103
HALT::dMinSecObl
double dMinSecObl
Halt integration at this obliquity.
Definition: eqtide.h:158
HALT::bPosDeDt
int bPosDeDt
Halt when de/dt > 0.
Definition: eqtide.h:161
SECONDARY::dDeDt
double dDeDt
Time derivative of eccentricity.
Definition: eqtide.h:125
PRIMARY::dSpinRate
double dSpinRate
Rotation Rate.
Definition: eqtide.h:85
SECONDARY
Definition: eqtide.h:101
SECONDARY::dTau
double dTau
Time lag.
Definition: eqtide.h:104
HALT::bSecLock
int bSecLock
Halt secondary becomes tide locked?
Definition: eqtide.h:163
HALT::bPriLock
int bPriLock
Halt primary becomes tide locked?
Definition: eqtide.h:162
PARAM::dMinValue
double dMinValue
Minimum value for e and obl.
Definition: eqtide.h:187
PARAM::bDoBackward
int bDoBackward
Perfrom backward integration?
Definition: eqtide.h:196
HALT::dMinSemi
double dMinSemi
Halt integration at this semi-major axis.
Definition: eqtide.h:156
PARAM::bLog
int bLog
Write log file?
Definition: eqtide.h:179
SECONDARY::dK2
double dK2
Secondary Love number.
Definition: eqtide.h:106
SECONDARY::dSemi
double dSemi
Semi-major axis.
Definition: eqtide.h:115
OPTIONS
Definition: eqtide.h:221
FILES::cIn
char * cIn
Input file.
Definition: eqtide.h:215
HALT::dMinEcc
double dMinEcc
Halt integration at this eccentricity.
Definition: eqtide.h:160
PARAM::iNumCols
int iNumCols
Number if columns in output file.
Definition: eqtide.h:185
IO::exit_input
int exit_input
Exit code for incompatible paramters.
Definition: eqtide.h:136
PARAM::iUnitTime
int iUnitTime
Time Unit for input/output.
Definition: eqtide.h:173
PARAM::dForwOutputTime
double dForwOutputTime
Output interval in forward integration (years)
Definition: eqtide.h:193
PARAM::iUnitAngle
int iUnitAngle
Angle Unit for input/output.
Definition: eqtide.h:172
PRIMARY::bForceEqSpin
int bForceEqSpin
Force spin rate to be equilibrium?
Definition: eqtide.h:89
PARAM::iIntegration
int iIntegration
Integration method.
Definition: eqtide.h:204
PARAM::iUnitLength
int iUnitLength
Length Unit for input/output.
Definition: eqtide.h:171
SECONDARY::dObliquity
double dObliquity
Obliquity.
Definition: eqtide.h:107
PRIMARY::dRadius
double dRadius
Radius of primary.
Definition: eqtide.h:80
PARAM::dForwStopTime
double dForwStopTime
Forward integration time (years)
Definition: eqtide.h:192
HALT::dMaxEcc
double dMaxEcc
Halt integration at this eccentricity.
Definition: eqtide.h:159
PRIMARY::iEpsilon
double * iEpsilon
Signs of phase lags.
Definition: eqtide.h:87
PRIMARY::iMassRad
int iMassRad
Mass-Radius relationship identifier.
Definition: eqtide.h:92
PARAM::iUnitMass
int iUnitMass
Mass Unit for input/output.
Definition: eqtide.h:170
PARAM::dl
double dl
Change in angular momentum: (l-l0)/l0.
Definition: eqtide.h:182
IO::exit_output
int exit_output
Exit code for error during output.
Definition: eqtide.h:139
PRIMARY::dObliquity
double dObliquity
Obliquity.
Definition: eqtide.h:84
SECONDARY::dEcc
double dEcc
Eccentricity.
Definition: eqtide.h:116
PARAM::bVarDt
int bVarDt
Use variable timestep?
Definition: eqtide.h:201
HALT::bHalt
int bHalt
Check for halts?
Definition: eqtide.h:153
IO::exit_io
int exit_io
Exit code for failure to open/close a file.
Definition: eqtide.h:135
SECONDARY::dMeanMotion
double dMeanMotion
Mean Motion.
Definition: eqtide.h:117
SECONDARY::dDomegaDt
double dDomegaDt
Time derivative of spin rate.
Definition: eqtide.h:122
IO::iDigits
int iDigits
Number of digits after the decimal.
Definition: eqtide.h:131
HALT::dMinPriObl
double dMinPriObl
Halt integration at this obliquity.
Definition: eqtide.h:157
PRIMARY::dMass
double dMass
Mass of primary.
Definition: eqtide.h:79
PARAM::fEqSpin
fdEqSpin fEqSpin
Function pointer to equilibrium spin rate.
Definition: eqtide.h:206
HALT::bMerge
int bMerge
Halt for merge?
Definition: eqtide.h:155
PARAM::dBackTimeStep
double dBackTimeStep
Time step for backward integration (years)
Definition: eqtide.h:197
PARAM::iTideModel
int iTideModel
Integer version of tide model.
Definition: eqtide.h:176
PRIMARY::dMaxLockDiff
double dMaxLockDiff
When to set spin rate to equilibrium.
Definition: eqtide.h:90
PARAM::dTimestepCoeff
double dTimestepCoeff
Variable Timestep coefficients.
Definition: eqtide.h:202
PRIMARY::dRG
double dRG
Radius of Gyration.
Definition: eqtide.h:86
PARAM::dForwTimeStep
double dForwTimeStep
Time step for forward integration (years)
Definition: eqtide.h:191
FILES
Definition: eqtide.h:213
PARAM::bDoForward
int bDoForward
Perform forward integration?
Definition: eqtide.h:190
PRIMARY::dTau
double dTau
Time lag.
Definition: eqtide.h:81
PRIMARY::dQ
double dQ
Tidal Q of primary.
Definition: eqtide.h:82
SECONDARY::dSpinRate
double dSpinRate
Rotation rate.
Definition: eqtide.h:108
FILES::cExe
char * cExe
Name of executable.
Definition: eqtide.h:214
PRIMARY
Definition: eqtide.h:78
HALT
Definition: eqtide.h:152
PARAM::l0
double l0
Initial angular momentum.
Definition: eqtide.h:181
SECONDARY::dMaxLockDiff
double dMaxLockDiff
When to set spin rate to equilibrium.
Definition: eqtide.h:113
PARAM::fDerivs
fvDerivs fDerivs
Function pointer to derivatives.
Definition: eqtide.h:205
Generated by   doxygen 1.8.9.1