# body¶

Contents

## body.c¶

Relationships between parameters associated with individual bodies.

This file contains subroutines that describe physical properties of any body. This includes conversions between the option parameter (a property that may be used at input) and the system parameter (the property in the BODY struct that is always up-to-date). If unsure between here and system.c, put here. Also includes mathemtatical relationships.

Author

Rory Barnes (RoryBarnes)

Date

May 7 2014

Functions

int fiSign(double dValue)

Calculate the sign of a number

Parameters
• dValue – The number whose sign is to be calculated

• iSign – The sign of the number

Returns

the sign (-1,0, or 1).

double fdDPerDt(double dRotRate, double dDrotrateDt)

Calcaulte the derivative of a period.

Parameters
• dRotRate – The rotational frequency

• dDrotrateDt – The time derivative of the rotational frequency

Returns

The time derivative of the rotational frequency

double fdTimescale(double dVar, double dDeriv)

Caclulate the characteristic timescale for a variable to change that is controlled by only 1 module.

Parameters
• dVar – The value of the variable

• dDeriv – The value of the variable’s time derivative

Returns

The timescale of the variable’s change: |x/(dx/dt)|. If the derivative is 0, return 0.

double fdTimescaleMulti(double dVar, double *dDeriv, int iNum)

Caclulate the characteristic timescale for a variable to change that is controlled by multiple processes.

Parameters
• dVar – The value of the variable

• dDeriv – Array of the values of the variable’s time derivatives

• iNum – The number of derivatives

• dTime – Dummy variable to keep track of the sum of the derivative

• iPert – Index of the perturbing process

Returns

The timescale of the variable’s change: |x/Sum(dx/dt)|

double fdFreqToPer(double dFreq)

Convert an angular frequency to a period

Parameters

dFreq – The frequency

Returns

The period

double fdPerToFreq(double dPeriod)

Convert a period to an angular frequency

Parameters

dPeriod – The period

Returns

The frequency

Calculate a body’s gravitational potential energy

Parameters
• dMass – The body’s mass

Returns

The body’s potential energy

Calculate a body’s rotational angular momentum

Parameters

• dMass – Body’s mass

• dOmega – Body’s rotational frequency

Returns

Body’s rotational angular momentum

Calculate a body’s rotational kinetic energy

Parameters
• dMass – Body’s mass

• dOmega – Body’s rotational frequency

Returns

Body’s rotational jinetic energy

Convert a rotational frequency to a rotational velocity

Parameters

• dFreq – Body’s rotational frequency

Returns

Body’s rotational velocity

Convert rotational velocity to rotational frequency

Parameters
• dRotVel – Body’s rotational velocity

Returns

Body’s rotational frequency

Calculate radius from density and mass

Parameters
• dDensity – Body’s bulk density

• dMass – Body’s mass

Returns

Calculate Mass from Radius and density

Parameters

• dDensity – Body’s bulk density

Returns

Body’s mass

Calculate rotational velocity from radius and rotation frequency

Parameters

• dRotRate – Body’s rotational frequency

Returns

Body’s rotational velocity

Calculate density for a uniform sphere

Parameters
• dMass – Body’s mass

Returns

Body’s bulk density

Stellar mass-radius relationship from New Light on Dark Stars, Table 4.1. See Barnes et al. (2013) Astrobiology 13:225-250.

Parameters

• x – Log base-10 of the stellar radius in solar units

• y – Log base-10 of the stellar mass in solar units

Returns

Stellar Mass

Stellar mass-radius relationship from New Light on Dark Stars, Table 4.1. See Barnes et al. (2013) Astrobiology 13:225-250.

Parameters
• dMass – Stellar mass

• x – Log base-10 of the stellar radius in solar units

• y – Log base-10 of the stellar mass in solar units

Returns

Stellar mass-radius relationship from Gorda, S. Yu. & Svechnikov, M. A. 1999, Astronomy Reports, 43, 521-525.

Parameters

dMass – Stellar Mass

Returns

Stellar mass-radius relationship from Gorda, S. Yu. & Svechnikov, M. A. 1999, Astronomy Reports, 43, 521-525. Reverse fit from Barnes et al. (2013) Astrobiology 13:225-250.

Parameters

• x – Log base-10 of the stellar radius in solar units

• y – Log base-10 of the stellar mass in solar units

Returns

Stellar mass

Stellar mass-radius relationship from Bayless, A.J. & Orosz, J.A. 2006, ApJ, 651, 1155-1165.

Parameters
• dMass – Stellar Mass

Returns

Stellar mass-radius relationship from Bayless, A.J. & Orosz, J.A. 2006, ApJ, 651, 1155-1165.

Parameters

• dMasss – Stellar masss is solar units

Returns

Stellar masss

Terrestrial planet mass-radius relationship from Sotin et al 2007, Icarus, 191, 337-351.

Parameters

dMass – Planetary mass

Returns

Terrestrial planet mass-radius relationship from Sotin et al 2007, Icarus, 191, 337-351.

Parameters

Returns

Planetary mass

Assign Radius based on mass and published relationship

Parameters
• dMass – Body’s mass

• iRelation – ID of mass-radius relationship to use

Returns

Body’s radius as provided by appropriate relationship

void BodyCopy(BODY *dest, BODY *src, EVOLVE *evolve)

Copy the body struct from src to dest

Parameters
• dest – Struct to receive the src

• src – Struct that contains original information

void CalcXYZobl(BODY *body, int iBody)

Calculate rotational variables from obliquity and precession angle

Parameters
• body – Body struct

• iBody – Index of the body struct for the body’s whose rotational spins is to be calculated.

double CalcDynEllipEq(BODY *body, int iBody)

Calculate equilibrium shape of planet using scaling laws and solar system values. If the value is less then Venus’, return Venus’.

Parameters
• bodyBODY struct

• iBody – Index of the body struct for the body’s whose equilibrium shape is to be calculated.

• J2Earth – Earth’s current oblateness

• J2Venus – Venus’ current oblateness

• CEarth – Earth’s current moment of inertia?

• nuEarth – ???

• EdEarth – ???

• dTmp – ???

• dDynEllip – Dynamical ellipticity

Returns

Dynamical elliptiticy

Lehmer+ (2017)’s model for the radius of a planet where it’s losing its atmopshere due to XUV radiation.

Parameters
• dMassEnv – Envelope’s mass

• dGravAccel – Body’s gravitational acceleration

• dScaleHeight – Atmospheric scale height

• dPresSurf – pressure at surface due to envelope

• dRadXUV – radius from center of planet where optical depth of XUV is unity

double fdLehmerPres(double dMassEnv, double dGravAccel, double dRadSurf)

Lehmer+ (2017)’s model for the pressure of a planet where it’s losing its atmopshere due to XUV radiation.

Parameters
• dMassEnv – Envelope’s mass

• dGravAccel – Body’s gravitational acceleration

• dPresXUV – Pressure at base of thermosphere

• dScaleHeight – Atmospheric scale height

• dPresSurf – pressure at surface due to envelope

double fdThermalTemp(BODY *body, int iBody)

Thermal temperature of an object heated by the radiation of its primary.

Parameters
• dFlux – Incident flux on the body

• dTemp – Thermal temperature

double fdImK2Total(BODY *body, int iBody)
double fdK2Man(BODY *body, int iBody)

Function compute upper mantle k2 Love number

Parameters
• body – Body struct

• iBody – Index of body

Returns

Upper mantle k2 Love number

double fdTidalQMan(BODY *body, int iBody)
double fdImK2ManThermint(BODY *body, int iBody)
double fdImK2Man(BODY *body, int iBody)

Function compute upper mantle imaginary component of k2 Love number

Parameters
• body – Body struct

• iBody – Index of body

Returns

Imaginary component of k2 Love number

double fdHflowSecMan(BODY *body, EVOLVE *evolve, int iBody)

Calculate total k_2 and Im(k_2) for a body. This value depends on which tidal model is called, and it the properties depend on the internal properties. void AssignTidalProperties(BODY *body,EVOLVE *evolve,int iBody) { if (body[iBody].bThermint) { The planet’s tidal response depends on mantle properties body[iBody].dK2=fdK2Man(body,iBody); body[iBody].dImK2=fdImK2Man(body,iBody); } else { We’re in a CPL/CTL type mode, and let’s start building the total K2 and ImK2 with the mantle. This should be sorted out in Verify. body[iBody].dK2 = body[iBody].dK2Man; body[iBody].dImK2 = body[iBody].dImK2Man; } Include tidal dissapation due to oceans if (body[iBody].bOceanTides) { weighted by the love number of each component body[iBody].dK2 += body[iBody].dK2Ocean; body[iBody].dImK2 += body[iBody].dImK2Ocean; } Include tidal dissapation due to envelope if (body[iBody].bEnvTides) { weighted by the love number of each component body[iBody].dK2 += body[iBody].dK2Env; body[iBody].dImK2 += body[iBody].dImK2Env; } Sanity checks: enforce upper bound if (body[iBody].dK2 > 1.5) { body[iBody].dK2 = 1.5; fprintf(stderr,”WARNING: body[%d].dK2 > 1.5 at time %.5e

years.\n”,iBody,evolve->dTime/YEARSEC); } } Function compute secular mantle heat flow: heat sinks - sources

Parameters
• body – Body struct

• iBody – Index of body

Returns

Heat flow of secular mantle cooling

int fiGetLowerBoundProximaCen(double dVal, const double *daArr, int iDim)

For use with `fdProximaCenStellar()` to interpolate stellar properties (temperature, radius, luminosity) from a grid.

Parameters
• dVal – Value of (temperature, radius, luminosity)

• daArr – Array of values from Yonsei-Yale tracks

• iDim – Length of array

• iIndex – Index of dArr correpoinding to the value

Returns

iIndex

double fdProximaCenBiLinear(int iALEN, double const data_lo[PROXIMACEN_ALEN], double const data_hi[PROXIMACEN_ALEN], int xi, int yi, double dx, double dy)

For use with `fdProximaCenStellar()` to interpolate stellar properties (temperature, radius, luminosity) from a grid. This function linearly interpolates over data, given indices of lower bounds on grid xi, yi and normalized distances to the interpolation point dx, dy.

What are these arguments?

double fdProximaCenInterpolate(int iALEN, int iMLEN, double const xarr[PROXIMACEN_ALEN], double const yarr[PROXIMACEN_MLEN], double const data_lo[PROXIMACEN_ALEN], double const data_hi[PROXIMACEN_ALEN], double A, double M, int *iError)

For use with `fdProximaCenStellar()` to interpolate stellar properties (temperature, radius, luminosity) from a grid.

What are these arguments?

double fdProximaCenStellar(int iParam, double A, double M, int *iError)

Computes the temperature, luminosity, or radius of Proxima Centauri by interpolating from a grid. The isochrones are from the Y^2 website

http://www.astro.yale.edu/demarque/fs255_grid.tar.gz

with [Fe/H] = +0.3 track and mixing length parameter alpha_MLT = 1.0. These were linearly interpolated between the 0.1 MSUN and 0.15 MSUN tracks to create luminosity, temperature, and radius interpolants, which are functions of a single variable (the age). The grid was then rectified and copied as constant arrays to `body.h`.

Note that in order to match the present-day luminosity (see below), I had to fudge a bit. I ended up scaling the entire luminosity track down by 15% to get it to match. I then re-computed the temperature to be consistent with the radius. Now all three quantities (L, R, and T) match within less than 1 sigma.

DATA FROM Boyajian+12; SECOND ROW FROM Demory+09 (direct measurements)

### 0.19 0.1410 ± 0.0070 0.00155 ± 0.00002 3054 ± 79¶

0.118 2.83E−04 0.00165 ± 0.00015 3098 ± 56 0.123 ± 0.006

What are these arguments?

int fiGetLowerBoundProximaCenB(double dVal, const double *daArr, int iDim)

For use with `fdProximaCenBRadius()` to interpolate the radius of Proxima Cen b from a grid, assuming it has a gaseous composition

Parameters
• dVal – Value of (temperature, radius, luminosity)

• daArr – Array of values from Yonsei-Yale tracks

• iDim – Length of array

• iIndex – Index of dArr correpoinding to the value

Returns

iIndex

double fdProximaCenBLinear(int xi, int yi, double dx, double dy)

For use with `fdProximaCenBRadius()` to interpolate the radius of Proxima Cen b from a grid, assuming it has a gaseous composition

What are the arguments?

double fdProximaCenBRadius(double C, double A, double M)

For use with `fdProximaCenBRadius()` to interpolate the radius of Proxima Cen b from a grid, assuming it has a gaseous composition

Here I’m assuming a mass of 1.27 MEARTH and a solid body radius of 1.074 REARTH. I’m using the Lopez+12 grids from Luger et al. (2015) and smoothing over sharp (presumably) numerical features.

What are the arguments?

double fdLopezRadius(double dMass, double dComp, double dFlux, double dAge, int iMetal)

Planet radius evolution from the Lopez et al. (2012) evolution grids.

What are the arguments?

double fdDotProduct(const int *x, const double *y)

Dot product of two vectors

Parameters
• x – First array

• y – Second array

• res – dot product

Returns

dot product of arrays x and y

void fvMatrixVectorMult(const int mat, const double *vec, double *result)

Matrix-vector multiplication

Parameters
• mat – Matrix

• vec – Vector

• result – Resultant vector

int fiGetLowerBound(double val, const double *arr, int dim)

Helper function for interpolating Baraffe grid

What are the arguments?

double fdBaraffeBiLinear(int iMLEN, int iALEN, double const data[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN], int xi, int yi, double dx, double dy)

Helper function for interpolating Baraffe grid

What are the arguments?

double fdBaraffeBiCubic(int iMLEN, int iALEN, double const data[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN], int xi, int yi, double dx, double dy)

Helper function for interpolating Baraffe grid

What are the arguments?

double fdBaraffeInterpolate(int iMLEN, int iALEN, double const xarr[STELLAR_BAR_MLEN], double const yarr[STELLAR_BAR_ALEN], double const data[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN], double M, double A, int iOrder, int *iError)

Helper function for interpolating Baraffe grid

What are the arguments?

double fdBaraffe(int iParam, double A, double M, int iOrder, int *iError)

Returns the stellar T, L, or R by interpolating over the Baraffe grid using either a bilinear (iOrder = 1) or a bicubic (iOrder = 3) interpolation.

What are the arguments?

void fdHabitableZoneKopparapu2013(BODY *body, int iNumBodies, double *daHZLimit)

Compute habitable zone limits from Kopparapu et al. (2013). Works with any number of stars.

Parameters
• body – An array of BODY structs

• iNumBodies – Total number of bodies in a system

• daHZLimit – A pointer to the array corresponding with the 6 limits

## body.h¶

Relationships between parameters associated with individual bodies.

Author

Rory Barnes (RoryBarnes)

Date

May 7 2014

Defines

REIDHAWLEY
GORDASVECH99
BAYLESSOROSZ06
SOTIN07
PROXIMACEN_T
PROXIMACEN_L
PROXIMACEN_R
PROXIMACEN_ERROR
PROXIMACEN_ALEN
PROXIMACEN_MLEN
PROXIMACEN_FUDGE
PROXCENBCOMPLEN
PROXCENBTIMELEN
MASSLEN
COMPLEN
FLUXLEN
METLEN
TIMELEN
STELLAR_T
STELLAR_L
STELLAR_R
STELLAR_RG
STELLAR_ERR_LINEAR
STELLAR_ERR_NONE
STELLAR_ERR_OUTOFBOUNDS_LO
STELLAR_ERR_OUTOFBOUNDS_HI
STELLAR_ERR_ISNAN
STELLAR_ERR_FILE
STELLAR_BAR_MLEN
STELLAR_BAR_ALEN

Variables

static double const PROXIMACEN_MARR[PROXIMACEN_MLEN] = {0.1, 0.15}
static double const PROXIMACEN_AARR[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGL_LO[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGR_LO[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGT_LO[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGL_HI[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGR_HI[PROXIMACEN_ALEN]
static double const PROXIMACEN_LOGT_HI[PROXIMACEN_ALEN]
static double const STELLAR_BAR_MARR[STELLAR_BAR_MLEN] = {0.0698, 0.07, 0.072, 0.075, 0.08, 0.09, 0.1, 0.11, 0.13, 0.15, 0.17, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.41}
static double const STELLAR_BAR_AARR[STELLAR_BAR_ALEN]
static const int STELLAR_BICUBIC_MATRIX = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, {-3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0}, {9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1}, {-6, 6, 6, -6, -3, -3, 3, 3, -4, 4, -2, 2, -2, -2, -1, -1}, {2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0}, {-6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1}, {4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1}}
static double const DATA_LOGL[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN]