ppo.h

Main header for the photodynamical routines.

Defines

MDFAST

Use the Murray & Dermott fast Kepler solver

NEWTON

Use the standard Newton Kepler solver

QGSL

Use the GSL complex polynomial solver

MAP_NONE

No map defined

MAP_RADIAL_DEFAULT

Default radially-symmetric LimbDarkenedMap

MAP_RADIAL_CUSTOM

Custom compiled Python radially-symmetric map

MAP_ELLIPTICAL_DEFAULT

Default elliptically-symmetric RadiativeEquilibriumMap

MAP_ELLIPTICAL_CUSTOM

Custom compiled Python elliptically-symmetric map

ERR_NONE

No error occurred

ERR_NOT_IMPLEMENTED

Function/option not yet implemented

ERR_KEPLER

Error in the Kepler solver; probably didn’t converge

ERR_INPUT

Bad input value

ERR_TOO_FEW_OCCS

No or too few occultations detected over the specified time interval

ERR_OOB

Warning: out of bounds

PI

Plain old pi

BIGG

Gravitational constant (mks)

DAYSEC

Number of seconds in one day

KBOLTZ

Boltzmann constant in W s / K

SBOLTZ

Stefan-Boltzmann constant in W / m^2 / K^4

SEARTH

Solar constant in W / m^2

HPLANCK

Planck constant in m^2 kg / s

CLIGHT

Speed of light in m / s

REARTH

Radius of the Earth in m

MICRON

Meters in 1 micron

PARSEC

Meters in 1 parsec

MEARTH

Mass of Earth in kg

GEARTH

Graviational constant in Earth units

MAXIM

Maximum magnitude of the imaginary component for root to be treated as real

SMALL

Tolerance in the geometry routines

TINY

Tolerance in the geometry routines

MINCRESCENT

Numerical issues pop up when the dayside crescent is too thin, so for phase angles smaller than this we use more zenith slices

CRESCENTNZ

Number of zenith slices in the tiny crescent limit

Typedefs

typedef double (*RADIANCEMAP)(double, double)

A radiance map function of the wavelength and the zenith angle.

typedef double (*CURVE)(double, ELLIPSE *)

A generic ellipse function.

typedef double (*INTEGRAL)(double, double, ELLIPSE *, int *)

A generic ellipse integral.

Functions

void GetAngles(double x, double y, double z, double vx, double vy, double vz, double Lambda, double Phi, double *theta, double *gamma)

Computes the phase angle theta and the rotation angle gamma at a given point in a planet’s orbit.

Parameters

• x: The x coordinate of the planet’s position on the sky
• y: The y coordinate of the planet’s position on the sky
• z: The z coordinate of the planet’s position on the sky
• vx: The x coordinate of the planet’s velocity on the sky
• vy: The y coordinate of the planet’s velocity on the sky
• vz: The z coordinate of the planet’s velocity on the sky
• Lambda: The latitudinal hotspot offset in radians (north positive)
• Phi: The longitudinal hotspot offset in radians (east positive)
• theta: The eyeball phase angle
• gamma: The eyeball rotation angle

double Blackbody(double lambda, double T)

Computes the blackbody intensity (W / m^2 / m / sr).

Return

The intensity given by Planck’s Law

Parameters

• lambda: The wavelength in m
• T: The effective temperature in K

int NextOccultation(int nt, double time[nt], int np, BODY **body, SETTINGS settings, int occulted, int noccultors, int occultors[noccultors], int noccultations, double occultation_times[noccultations], int occultation_inds[noccultations], double occultation_durs[noccultations])

Computes the time(s) of the next occultation(s) of the body of index occulted by the body(ies) of index(ices) occultors.

Return

The error code

Parameters

• nt: The size of the time array
• time: The array of times at which to evaluate the orbits
• np: The number of particles (bodies) in the system
• body: An array of BODY pointers corresponding to all the bodies in the system
• settings: An instance of the SETTINGS class containing all settings
• occulted: The index of the body being occulted
• noccultors: Size of occultor array
• occultors: The indices of the occultors of occulted
• noccultations: The number of occultations to look for
• occultation_times: The times of each of the occultations, computed by this function
• occultation_inds: The indices of the occultors corresponding to each of the occultations identified by this function
• occultation_durs: The duration of each of the occultations identified by this function

int NBody(int np, BODY **body, SETTINGS settings, int halt_on_occultation, int occulted, int noccultors, int occultors[noccultors], int noccultations, double occultation_times[noccultations], int occultation_inds[noccultations], double occultation_durs[noccultations])

Compute the instantaneous x, y, and z positions of all planets using the REBOUND N-Body code.

Return

The error code

Parameters

• np: The number of particles (bodies) in the system
• body: An array of BODY pointers to each of the bodies in the system
• settings: An instance of the SETTINGS class
• halt_on_occultation: Halt when a certain number of occultations of a given body have occurred?
• occulted: The index of the body being occulted (only used if halt_on_occultation = 1)
• noccultors: Size of occultor array (only used if halt_on_occultation = 1)
• occultors: The indices of the occultors of occulted (only used if halt_on_occultation = 1)
• noccultations: The number of occultations to look for (only used if halt_on_occultation = 1)
• occultation_times: The times of each of the occultations, computed by this function (only used if halt_on_occultation = 1)
• occultation_inds: The indices of the occultors corresponding to each of the occultations identified by this function (only used if halt_on_occultation = 1)
• occultation_durs: The duration of each of the occultations identified by this function (only used if halt_on_occultation = 1)

int Kepler(int np, BODY **body, SETTINGS settings)

Compute the instantaneous x, y, and z positions of all planets with a simple multi-Keplerian solver.

Return

The error code

Parameters

• np: The number of particles (bodies) in the system
• body: An array of BODY pointers to each of the bodies in the system
• settings: An instance of the SETTINGS class

void OccultedFlux(double r, int no, double x0[no], double y0[no], double ro[no], double theta, double tnight, double teff, double distance, double mintheta, int maxvertices, int maxfunctions, int adaptive, int circleopt, int batmanopt, int quarticsolver, int nu, int nz, int nw, double u[nu *nw], double lambda[nw], double flux[nw], int maptype, RADIANCEMAP radiancemap, int quiet, int *iErr)

Computes the flux of an eyeball planet occulted by one or more bodies over a grid of wavelengths.

Parameters

• r: The radius of the occulted planet
• no: The number of occulting bodies
• x0: The x coordinate of the center of each occulting body
• y0: The y coordinate of the center of each occulting body
• ro: The radius of each occulting body
• theta: The eyeball phase angle
• tnight: The planet’s night side temperature in K (airless limit only)
• teff: The planet’s effective temperature in K (blackbody limit only)
• distance: The distance to the system in parsecs
• mintheta: The minimum absolute value of the phase angle, below which it is assumed constant to prevent numerical errors
• maxvertices: Maximum number of vertices in the problem
• maxfunctions: Maximum number of functions in the problem
• adaptive: Adaptive zenith angle grid? Limb-darkened limit only
• circleopt: Treat ellipses as circles for limb-darkened bodies? There’s no reason not to do this! Option left in mainly for testing purposes.
• batmanopt: Use the BATMAN algorithm to speed up limb-darkened occultations?
• quarticsolver: Which quartic solver to use?
• nu: The number of limb darkening coefficients
• nz: The number of zenith angle slices
• nw: The size of the wavelength grid
• u: The limb darkening coefficient grid: nu coefficients at each wavelength
• lambda: The wavelength grid in m
• flux: The wavelength-dependent light curve computed by this function
• maptype: The code for the radiance map
• radiancemap: A pointer to the RADIANCEMAP function
• quiet: Suppress output?
• iErr: Flag set if an error occurs

void UnoccultedFlux(double r, double theta, double tnight, double teff, double distance, double mintheta, int maxvertices, int maxfunctions, int adaptive, int circleopt, int batmanopt, int quarticsolver, int nu, int nz, int nw, double u[nu *nw], double lambda[nw], double flux[nw], int maptype, RADIANCEMAP radiancemap, int quiet, int *iErr)

Computes the total flux of an eyeball over a grid of wavelengths. This routine (hackishly) calls OccultedFlux with an imaginary occultor covering the entire disk of the body.

Parameters

• r: The radius of the occulted planet
• theta: The eyeball phase angle
• tnight: The planet’s night side temperature in K (airless limit only)
• teff: The planet’s effective temperature in K (blackbody limit only)
• distance: The distance to the system in parsecs
• mintheta: The minimum absolute value of the phase angle, below which it is assumed constant to prevent numerical errors
• maxvertices: Maximum number of vertices in the problem
• maxfunctions: Maximum number of functions in the problem
• adaptive: Adaptive zenith angle grid? Limb-darkened limit only
• circleopt: Treat ellipses as circles for limb-darkened bodies? There’s no reason not to do this! Option left in mainly for testing purposes.
• batmanopt: Use the BATMAN algorithm to speed up limb-darkened occultations?
• quarticsolver: Which quartic solver to use?
• nu: The number of limb darkening coefficients
• nz: The number of zenith angle slices
• nw: The size of the wavelength grid
• u: The limb darkening coefficient grid: nu coefficients at each wavelength
• lambda: The wavelength grid in m
• flux: The wavelength-dependent light curve computed by this function
• maptype: The code for the radiance map
• radiancemap: A pointer to the RADIANCEMAP function
• quiet: Suppress output?
• iErr: Flag set if an error occurs

int Orbits(int nt, double time[nt], int np, BODY **body, SETTINGS settings)

Evolves the orbital positions of all bodies forward in time with either a Keplerian or an N-body solver.

Return

The error code

Parameters

• nt: The size of the time array
• time: The array of times at which to evaluate the orbits
• np: The number of particles (bodies) in the system
• body: An array of BODY pointers corresponding to all the bodies in the system
• settings: An instance of the SETTINGS class containing all settings

int Flux(int nt, double time[nt], int nw, double wavelength[nw], double continuum[nt *nw], int np, BODY **body, SETTINGS settings)

Computes the full light curve for all bodies in the system.

Return

The error code

Parameters

• nt: The size of the time array
• time: The array of times at which to evaluate the orbits
• nw: The size of the wavelength grid
• wavelength: The wavelength grid in microns
• continuum: The continuum flux (i.e., the total flux without occultations) of the system on a time/wavelength grid
• np: The number of particles (bodies) in the system
• body: An array of BODY pointers corresponding to all the bodies in the system
• settings: An instance of the SETTINGS class containing all settings

struct BODY

#include <ppo.h>

Struct containing all the information pertaining to a body (star, planet, or moon) in the system.

Public Members

double m

Body mass in Earth masses

double per

Orbital period in days

double inc

Orbital inclination in radians

double ecc

Orbital eccentricity

double w

Longitude of pericenter in radians

double Omega

Longitude of ascending node in radians

double a

Orbital semi-major axis in Earth radii

double tperi0

Time of pericenter passage in days

double r

Body radius in Earth radii

double albedo

Body albedo

double teff

Effective temperature in K

double tnight

Night side temperature in K

int phasecurve

Compute the phase curve for this body?

double Lambda

Longitudinal hotspot offset in radians

double Phi

Latitudinal hotspot offset in radians

int host

The index of this body’s host

int cartesian

Initialize body based on its Cartesian coords?

int nu

Number of limb darkening coefficients (per wavelength)

int nz

Number of zenith angle slices

int nt

Size of time array

int nw

Size of wavelength array

double *u

Limb darkening coefficient grid

double *time

Time array

double *wavelength

Wavelength array

double *x

The Cartesian x position on the sky (right positive)

double *y

The Cartesian y position on the sky (up positive)

double *z

The Cartesian z position on the sky (into sky positive)

double *vx

The Cartesian x velocity on the sky (right positive)

double *vy

The Cartesian y velocity on the sky (up positive)

double *vz

The Cartesian z velocity on the sky (into sky positive)

int *occultor

The array of occultor bit flags

double *flux

The grid of observed flux from this body in time/wavelength

double *total_flux

The total unocculted flux of this body at full phase

int maptype

Which radiance map to use?

RADIANCEMAP radiancemap

A function that returns a radiance map for the body’s surface

struct SETTINGS

#include <ppo.h>

Struct containing all the settings for computing orbits and light curves for the system.

Public Members

int nbody

Use N-Body solver?

int integrator

Which N-Body solver to use

double keptol

Kepler solver tolerance

int maxkepiter

Maximum iterations in Kepler solver

int kepsolver

Which Kepler solver to use

double timestep

N-Body timestep in days

int adaptive

Adaptive zenith angle grid for limb-darkened bodies?

int circleopt

Treat zenith angle slices as circles in the limb-darkened limit? No reason to set this to 0

int batmanopt

Use BATMAN algorithm to speed up limb-darkened occultation light curves?

int quarticsolver

Which quartic solver to use?

int quiet

Suppress output?

double mintheta

Minimum absolute value of the phase angle below which it is assumed to be constant to prevent numerical errors

int maxvertices

Maximum number of vertices (for memory allocation)

int maxfunctions

Maximum number of functions (for memory allocation)

double distance

Distance to the system in parsecs

int nstars

Number of stars in the system

struct ELLIPSE

#include <ppo.h>

A generic ellipse class.

Public Members

int circle

Is this a circle?

double r

Radius (if it’s a circle)

double a

Semi-major axis (if it’s an ellipse)

double b

Semi-minor axis (if it’s an ellipse)

double x0

x coordinate of the center of the ellipse

double y0

y coordinate of the center of the ellipse

double xmin

Leftmost point of the ellipse

double xmax

Rightmost point of the ellipse

struct FUNCTION

#include <ppo.h>

A container for the properties, functions, and integrals of an ellipse.

Public Members

double y

The value of the function at a given point

ELLIPSE *ellipse

The ellipse class corresponding to this function

CURVE curve

The curve function

INTEGRAL integral

The integral function