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

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]