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

double fdBodyPotEnergy(double dMass, double dRadius)

Calculate a body’s gravitational potential energy

Parameters:

dMass – The body’s mass
dRadius – The body’s radius

Returns:

The body’s potential energy

double fdRotAngMom(double dRadGyra, double dMass, double dRad, double dOmega)

Calculate a body’s rotational angular momentum

Parameters:

dRadGyra – Body’s radius of gyration
dMass – Body’s mass
dRad – Body’s radius
dOmega – Body’s rotational frequency

Returns:

Body’s rotational angular momentum

double fdRotKinEnergy(double dMass, double dRadius, double dRadGyra, double dOmega)

Calculate a body’s rotational kinetic energy

Parameters:

dMass – Body’s mass
dRadius – Body’s radius
dRadyGyra – Body’s radius of gyration
dOmega – Body’s rotational frequency

Returns:

Body’s rotational jinetic energy

double fdRadiusFreqToRotVel(double dRadius, double dFreq)

Convert a rotational frequency to a rotational velocity

Parameters:

dRadius – Body’s radius
dFreq – Body’s rotational frequency

Returns:

Body’s rotational velocity

double fdRadiusRotVelToFreq(double dRotVel, double dRadius)

Convert rotational velocity to rotational frequency

Parameters:

dRotVel – Body’s rotational velocity
dRadius – Body’s radius

Returns:

Body’s rotational frequency

double fdDensityMassToRadius(double dDensity, double dMass)

Calculate radius from density and mass

Parameters:

dDensity – Body’s bulk density
dMass – Body’s mass

Returns:

Body’s Radius

double fdMassFromRadiusDensity(double dRadius, double dDensity)

Calculate Mass from Radius and density

Parameters:

dRadius – Body’s radius
dDensity – Body’s bulk density

Returns:

Body’s mass

double fdRotVel(double dRadius, double dRotRate)

Calculate rotational velocity from radius and rotation frequency

Parameters:

dRadius – Body’s radius
dRotRate – Body’s rotational frequency

Returns:

Body’s rotational velocity

double fdSphereDensity(double dMass, double dRadius)

Calculate density for a uniform sphere

Parameters:

dMass – Body’s mass
Body's – radius

Returns:

Body’s bulk density

double fdRadToMass_ReidHawley(double dRadius)

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

Parameters:

dRadius – Stellar radius
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

double fdMassToRad_ReidHawley(double dMass)

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 Radius

double fdMassToRad_GordaSvech99(double dMass)

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

Parameters:: dMass – Stellar Mass
Returns:: Stellar radius

double fdRadToMass_GordaSvech99(double dRadius)

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:

dRadius – Stellar Radius
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

double fdMassToRad_BaylessOrosz06(double dMass)

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

Parameters:

dMass – Stellar Mass
dRadius – Stellar radius is solar units

Returns:

Stellar radius

double fdRadToMass_BaylessOrosz06(double dRadius)

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

Parameters:

dRadius – Stellar Radius
dMasss – Stellar masss is solar units

Returns:

Stellar masss

double fdMassToRad_Sotin07(double dMass)

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

Parameters:: dMass – Planetary mass
Returns:: Planetary radius

double fdMassToRad_LehmerCatling17(double dMass)

double fdRadToMass_Sotin07(double dRadius)

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

Parameters:: dRadius – Planetary radius
Returns:: Planetary mass

double fdMassToRad(double dMass, int iRelation)

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

double fdRadToMass(double dMass, int iRelation)

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:

body – BODY 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

double fdLehmerRadius(BODY *body, int iBody)

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
dRadSurf – Surface Radius of rocky core
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
dRadSurf – Surface Radius of rocky core
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)

Fe/H Radius Luminosity Teff Mass Lx/Lbol

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[16][16], 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

double fdEffectiveTemperature(BODY *body, int iBody)

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_ERR_BADORDER

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[16][16] = {{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]

static double const DATA_RADIUS[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN]

static double const DATA_LOGT[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN]

static double const DATA_RG[STELLAR_BAR_MLEN][STELLAR_BAR_ALEN]

body

body.c

Fe/H Radius Luminosity Teff Mass Lx/Lbol

0.19 0.1410 ± 0.0070 0.00155 ± 0.00002 3054 ± 79

body.h

body.c 

body.h 