binary

Contents

binary.c 

Subroutines that control the integration of the circumbinary planet orbital dynamics module.

Author

David Fleming (dflemin3)

Date

Jan 12 2016

Description

Module to model circumbinary planet dynamics.

(circumbinary planet(s))

body approximation and hence the CBPs are not allowed to graviationally interact.

Note

body 0 = primary star, body 1 = secondary star, body 2+ = CBP

Note

The :cite:`Leung2013` theory ONLY applies to the restricted 3

Functions

void BodyCopyBinary(BODY *dest, BODY *src, int foo, int iNumBodies, int iBody): Copy body properties from src to dest for cbp

void InitializeBodyBinary(BODY *body, CONTROL *control, UPDATE *update, int iBody, int iModule): Only use this function for malloc’ing stuff Since nothing has to be malloc’ed for binary, do nothing

void InitializeUpdateTmpBodyBinary(BODY *body, CONTROL *control, UPDATE *update, int iBody): No need to allocate anything

void ReadFreeEcc(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): Free eccentricity, can’t be in primary file Note: Do error checking for negative values

void ReadLL13N0(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): Lee + Leung 2013 Mean Motion N0 This parameter cannot exist in primary file

void ReadLL13K0(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): Lee + Leung 2013 radial epicyclic frequency This parameter cannot exist in primary file

void ReadLL13V0(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): Lee + Leung 2013 vertical epicyclic frequency This parameter cannot exist in primary file

void ReadFreeInc(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in the primary file Free inclination goes from 0 to 180 degrees

void ReadLL13PhiAB(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in the primary file PhiAB goes from [0,360) if in degrees

void ReadCBPM0(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in the primary file CBPM0 goes from [0,360) if in degrees

void ReadCBPZeta(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in the primary file. dCBPZeta goes from [0,360) if in degrees

void ReadCBPPsi(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in the primary file. dCBPPsi goes from [0,360) if in degrees

void ReadHaltHolmanUnstable(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in primary file

void ReadHaltRocheLobe(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile): This parameter cannot exist in primary file

void InitializeOptionsBinary(OPTIONS *options, fnReadOption fnRead[]): Initialization Options for BINARY

void ReadOptionsBinary(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, fnReadOption fnRead[], int iBody): Read all BINARY input options.

void VerifyCBPR(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void VerifyCBPZ(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void VerifyCBPPhi(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void VerifyCBPRDot(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void VerifyCBPZDot(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void VerifyCBPPhiDot(BODY *body, OPTIONS *options, UPDATE *update, double dAge, int iBody)

void fnPropsAuxBinary(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update, int iBody)

void fnForceBehaviorBinary(BODY *body, MODULE *module, EVOLVE *evolve, IO *io, SYSTEM *system, UPDATE *update, fnUpdateVariable ***fnUpdate, int iBody, int iModule)

void AssignBinaryDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, fnUpdateVariable ***fnUpdate, int iBody)

void NullBinaryDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, fnUpdateVariable ***fnUpdate, int iBody)

void VerifyBinary(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, OUTPUT *output, SYSTEM *system, UPDATE *update, int iBody, int iModule)

void InitializeUpdateBinary(BODY *body, UPDATE *update, int iBody)

void FinalizeUpdateCBPRBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

void FinalizeUpdateCBPZBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

void FinalizeUpdateCBPPhiBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

void FinalizeUpdateCBPRDotBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

void FinalizeUpdateCBPZDotBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

void FinalizeUpdateCBPPhiDotBinary(BODY *body, UPDATE *update, int *iEqn, int iVar, int iBody, int iFoo)

int fbHaltHolmanUnstable(BODY *body, EVOLVE *evolve, HALT *halt, IO *io, UPDATE *update, fnUpdateVariable ***fnUpdate, int iBody): If the CBP’s dSemi is less than the Holman stability limit, it’s unstable and integration ends

int fbHaltRocheLobe(BODY *body, EVOLVE *evolve, HALT *halt, IO *io, UPDATE *update, fnUpdateVariable ***fnUpdate, int iBody): If the secondary enters the roche lobe of the primary, HALT!

void CountHaltsBinary(HALT *halt, int *iHalt)

void VerifyHaltBinary(BODY *body, CONTROL *control, OPTIONS *options, int iBody, int *iHalt)

void WriteFreeEccBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteFreeIncBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteBinPriRBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[]): Write the binary primary star radial position

void WriteBinSecRBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[]): Write the binary secondary star radial position

void WriteBinPriPhiBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteBinSecPhiBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPPhiBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPPhiDotBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteLL13N0Binary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteLL13K0Binary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteLL13V0Binary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPRBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[]): Write the circumbinary planet orbital radius (CBPR)

void WriteCBPR0Binary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[]): Write the circumbinary planet guiding radius (CBPR0)

void WriteCBPZBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPRDotBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPZDotBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[])

void WriteCBPInsol(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char cUnit[]): Write Earth-normalized insolation received by CBP averaged over 1 binary orbit assuming P_bin << P_cbp

void InitializeOutputBinary(OUTPUT *output, fnWriteOutput fnWrite[])

void LogOptionsBinary(CONTROL *control, FILE *fp)

void LogBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UPDATE *update, fnWriteOutput fnWrite[], FILE *fp): Anything here? Nope.

void LogBodyBinary(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UPDATE *update, fnWriteOutput fnWrite[], FILE *fp, int iBody)

void AddModuleBinary(CONTROL *control, MODULE *module, int iBody, int iModule)

double fndDot(double *a, double *b): 3D dot product

int fniDelta(int i, int j): Kronecker Delta

void fnvCylToCartPos(double *daCylPos, double *dCartPos): Convert from cylindrical position coords to cartesian (3 dimensional)

void fnvCylToCartVel(double *daCylPos, double *daCylVel, double *dCartVel): Convert from cylindrical velocity coords to cartesian (3 dimensional)

void fnvSpecificAngMom(double *r, double *v, double *h): Calculate specific angular momentum h = r x v in cartesian coords

double fndSpecificOrbEng(BODY *body, int iBody): Calculate specific orbital energy eps = v^2/2 - mu/|r|

void fnvAssignOrbitalElements(BODY *body, int iBody): Compute and assign CBP’s orbital elements

double fndComputeSemi(BODY *body, int iBody): Compute CBP’s semi-major axis a = -mu/(2*eps)

double fndComputeEcc(BODY *body, int iBody): Compute CBP’s orbital eccentricity e = sqrt(1 + 2*eps*h*h/(mu*mu))

double fndComputeInc(BODY *body, int iBody): Compute CBP’s inclincation i = arccos(h_z/|h|)

double fndComputeLongA(BODY *body, int iBody): Compute a CBP’s longitude of ascending node See: https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node

void fnvComputeEccVector(BODY *body, double *evec, int iBody): Compute a CBP’s eccentricity vector

double fndComputeArgPeri(BODY *body, int iBody): Compute a CBP’s argument of pericenter

double fndMeanToEccentric(double M, double e)

Solves kepler’s equation via Newton’s method

Parameters:

M – double, Mean anomaly
e, eccentricity –

double fndEccToTrue(double E, double e): Convert eccentric anomaly to true anomaly

double fndHolmanStability(BODY *body): Compute the dynamical stability limit, a_crit, first computed by Holman and Wiegert 1999 that depends on binary dSemi, dEcc If CBP.dSemi < a_crit, planet is unstable and system should halt

double fndRocheLobe(BODY *body): Compute the roche lobe of the primary according to Egglton’s formula given in https://en.wikipedia.org/wiki/Roche_lobe

double fndBinaryMeanAnomaly(double dMeanMotion, double dTime, double dPhi): Compute the binary mean anomaly M = n*t + phi where phi is an arbitrary offset Note: When used with the binary, dPhi == dLL13PhiAB

double fndMeanMotionBinary(BODY *body, int iBody): LL13 N0

double fndEpiFreqK(BODY *body, int iBody): LL13 K0

double fndEpiFreqV(BODY *body, int iBody): LL13 V0

double fndPhi0(double dTime, double dMeanMotion, double dVarPhi): Circular (azimuthal) motion of the guiding center for a cbp: phi0 dVarPhi here is equal to dCBPM0, the initial CBP mean anomaly LL13 Eqn 20

double fndPot0(int j, int k, double R, BODY *body): LL13 Eqn 15: Binary potential component 0

double fndPot0dR(int j, int k, double R, BODY *body): LL13 Eqn 15: d/dR of binary potential component 0

double fndPot1(int j, int k, double R, BODY *body): LL13 eqn 16: Binary potential component 1

double fndPot1dR(int j, int k, double R, BODY *body): LL13 eqn 16: d/dR of binary potential component 1

double fndn(double R, BODY *body): LL13 function defined in paragraph in between eqn 24 and 25 which is like the mean motion at R

double fndC0(BODY *body, int iBody): LL13 C0 as defined in eqn 28

double fndC0k(int k, BODY *body, int iBody): LL13 C^0_k as defined by eqn 29

double fndCPk(int k, BODY *body, int iBody): LL13 C^+_k as defined by eqn 30a (the + term)

double fndCMk(int k, BODY *body, int iBody): LL13 C^-_k as defined by eqn 30b (the - term)

double fndD0(BODY *body, int iBody): LL13 D0 as defined by eqn 32

double fndDk0(int k, BODY *body, int iBody): LL13 Dk0 as defined by eqn 33

double fndDPk(int k, BODY *body, int iBody): LL13 D+_k as defined by eqn 34a (the + term)

double fndDMk(int k, BODY *body, int iBody): LL13 D-_k as defined by eqn 34b (the - term)

double fndCBPRBinary(BODY *body, SYSTEM *system, int *iaBody): Computes the CBP orbital radius

double fndCBPPhiBinary(BODY *body, SYSTEM *system, int *iaBody)

double fndCBPZBinary(BODY *body, SYSTEM *system, int *iaBody)

double fndCBPRDotBinary(BODY *body, SYSTEM *system, int *iaBody)

double fndCBPPhiDotBinary(BODY *body, SYSTEM *system, int *iaBody)

double fndCBPZDotBinary(BODY *body, SYSTEM *system, int *iaBody)

double fndFluxApproxBinary(BODY *body, int iBody): Compute the approximate flux in the limit P_bin << P_cbp received by the CBP from the 2 stars averaged over 1 binary orbit Assumes binary orb elements don’t vary much per orbit (safe assumption)

double fndFluxExactBinary(BODY *body, int iBody, double L0, double L1): Compute the exact flux (as close to exact as you want) received by the CBP from the 2 stars averaged over 1 CBP orbit Assumes binary orb elements don’t vary much over 1 CBP orbit Also assumes that 1 CBP orbit is approximately Keplerian

double fndApproxEqTemp(BODY *body, int iBody, double dAlbedo): Compute the approximate equlibrirum temperature for a circumbinary planet averaged over the binary orbit

double fndApproxInsol(BODY *body, int iBody): Compute the approximate insolation for a circumbinary planet averaged over the binary orbit relative to Earth’s where F_Earth = 1361 W/m^2

void fnvBinaryDebug(BODY *body): Dumps out a bunch of values to see if they agree with LL13

binary.h 

Subroutines that control the integration of the circumbinary planet orbital dynamics module.

Author: David Fleming (dflemin3)
Date: Jan 12 2016

Defines

K_MAX

FLUX_INT_MAX

KEQNTOL

MAX_KEPLER_ITERS

FLUX_EARTH

OPTSTARTBINARY

OPTENDBINARY

OPT_FREEECC

OPT_FREEINC

OPT_LL13N0

OPT_LL13K0

OPT_LL13V0

OPT_LL13PHIAB

OPT_CBPM0

OPT_CBPZETA

OPT_CBPPSI

OPT_HALTHOLMAN

OPT_HALTROCHELOBE

OUTSTARTBINARY

OUTENDBINARY

OUT_FREEECC

OUT_FREEINC

OUT_LL13N0

OUT_LL13K0

OUT_LL13V0

OUT_CYLPOS

OUT_CBPR

OUT_CBPZ

OUT_CBPPHI

OUT_CBPRDOT

OUT_CBPZDOT

OUT_CBPPHIDOT

OUT_CBPR0

OUT_CBPINSOL

OUT_BINPRIR

OUT_BINPRIPHI

OUT_BINSECR

OUT_BINSECPHI