eyeball.c

Occultation light curve routines.

Functions

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

void BatmanFlux(double r, double x0, double y0, double ro, double teff, double distance, int nu, int nz, int nw, double u[nu *nw], double lambda[nw], double flux[nw])

An implementation of the BATMAN algorithm (Kreidberg 2016) to compute occultations of radially-symmetric bodies. This applies to transits across stars and planet-planet occultations in the limb-darkened limit or of planets at full or new phase.

Parameters
  • r: The occulted body’s radius in Earth radii
  • x0: The x coordinate of the occultor relative to the occulted body in Earth radii
  • y0: The y coordinate of the occultor relative to the occulted body in Earth radii
  • ro: The radius of the occultor in Earth radii
  • teff: The effective temperature of the occulted body
  • distance: The distance to the system in parsecs
  • 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

int dblcomp(const void *a, const void *b)

Compares two doubles. Used with qsort.

Return
0 if a == b, 1 if a < b, -1 otherwise
Parameters
  • a: A double
  • b: A double

int funcomp(const void *a, const void *b)

Compares the y values of two FUNCTIONs. Used with qsort.

Return
0 if a.y == b.y, 1 if a.y < b.y, -1 otherwise
Parameters

double RadiativeEquilibriumMap(double lambda, double za, double teff, double tnight)

Computes the flux at a given zenith angle on an eyeball planet.

Return
The emitted spectral flux density from that point on the surface
Parameters
  • lambda: The wavelength in m
  • za: The zenith angle in radians, measured from the hotspot
  • teff: The body’s equilibrium temperature
  • tnight: The night side temperature in K

double LimbDarkenedMap(int j, double za, double teff, int nu, int nw, double u[nu *nw], double B0)

Computes the flux at a given zenith angle on a limb-darkened planet.

Return
The emitted spectral flux density from that point on the surface
Parameters
  • j: The index of the wavelength grid
  • za: The zenith angle in radians, measured from the hotspot
  • teff: The body’s effective temperature
  • nu: The number of limb-darkening coefficients per wavelength bin
  • nw: The number of wavelength bins
  • u: The limb-darkening coefficient grid
  • B0: The normalization constant for the radiance at this wavelength

void GetRootsGSL(double a, double b, double xE, double yE, double xC, double yC, double r, double roots[4])

Computes the real roots of the circle-ellipse intersection quartic equation. These roots are the x coordinates of the intersection point(s) relative to the center of the circle.

Parameters
  • a: The semi-major axis of the ellipse, which is aligned with the y axis
  • b: The semi-minor axis of the ellipse, which is aligned with the x axis
  • xE: The x coordinate of the center of the ellipse
  • yE: The y coordinate of the center of the ellipse
  • xC: The x coordinate of the center of the circle
  • yC: The y coordinate of the center of the ellipse
  • r: The radius of the circle
  • roots: The x coordinates of the point(s) of intersection

double ZenithAngle(double x, double y, double r, double theta)

Computes the zenith angle of a point (x, y) on the projected disk of a planet.

Return
The zenith angle in radians
Parameters
  • x: The x coordinate of the point
  • y: The y coordinate of the point
  • r: The radius of the planet
  • theta: The phase angle of the eyeball

void SurfaceIntensity(double tnight, double teff, int maptype, RADIANCEMAP radiancemap, int nz, double zenithgrid[nz], int nw, double lambda[nw], int nu, double u[nu *nw], double B0[nw], double B[nz+1][nw])

Computes the blackbody intensity at the center of each zenith_angle slice evaluated at a given array of wavelengths.

Parameters
  • tnight: The night side temperature in K (eyeball limit)
  • teff: The effective temperature of the planet (blackbody limit)
  • maptype: The code for the radiance map
  • radiancemap: A pointer to the RADIANCEMAP function
  • nz: The number of zenith angle slices
  • zenithgrid: The zenith angle grid
  • nw: The number of wavelength points
  • lambda: The wavelength array
  • nu: The number of limb darkening coefficients
  • u: The limb darkening coefficient grid: nu coefficients at each wavelength
  • B0: The wavelength-dependent normalization for the limb-darkened case
  • B: The blackbody intensity grid

double curve(double x, FUNCTION function)

Evaluates a FUNCTION instance.

Return
The value of the function at x
Parameters
  • x: The point at which to evaluate the function
  • function: A FUNCTION instance

double integral(double x0, double x1, FUNCTION function, int *oob)

Evaluates the definite integral of a FUNCTION instance.

Return
The definite integral of the function from x0 to x1
Parameters
  • x0: The lower integration limit
  • x1: The upper integration limit
  • function: A FUNCTION instance
  • oob: Flag set to 1 if the integration limits are out of bounds

double fupper(double x, ELLIPSE *ellipse)

The upper curve of an ellipse.

Return
The value of the upper curve of the ellipse at x.
Parameters
  • x: The point at which to evaluate the curve
  • ellipse: An ELLIPSE instance

double flower(double x, ELLIPSE *ellipse)

The lower curve of an ellipse.

Return
The value of the lower curve of the ellipse at x.
Parameters
  • x: The point at which to evaluate the curve
  • ellipse: An ELLIPSE instance

double iupper(double xL, double xR, ELLIPSE *ellipse, int *oob)

The area under the upper curve of an ellipse.

Return
The value of the definite integral of the upper curve of the ellipse
Parameters
  • xL: The left (lower) integration limit
  • xR: The right (upper) integration limit
  • ellipse: An ELLIPSE instance
  • oob: Flag set to 1 if the integration limits are out of bounds

double ilower(double xL, double xR, ELLIPSE *ellipse, int *oob)

The area under the lower curve of an ellipse.

Return
The value of the definite integral of the lower curve of the ellipse
Parameters
  • xL: The left (lower) integration limit
  • xR: The right (upper) integration limit
  • ellipse: An ELLIPSE instance
  • oob: Flag set to 1 if the integration limits are out of bounds

void AddZenithAngleEllipse(double zenith_angle, double r, int no, double x0[no], double y0[no], double ro[no], double theta, int quarticsolver, int maxvertices, int maxfunctions, double *vertices, int *v, FUNCTION *functions, int *f)

Computes the full geometry for a single zenith angle ellipse. Computes the functional form of the ellipse, its integral, and all points of intersection with the limb of the occulted planet and the limb of the occultor. Adds the vertices and functions to the running lists.

Parameters
  • zenith_angle: The zenith angle in radians
  • 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
  • quarticsolver: Which quartic solver to use?
  • theta: The phase angle of the eyeball
  • maxvertices: Maximum number of vertices in the problem
  • maxfunctions: Maximum number of functions in the problem
  • vertices: The list of vertices
  • v: The number of vertices so far
  • functions: The list of functions
  • f: The number of functions so far

void AddZenithAngleCircle(double zenith_angle, double r, int no, double x0[no], double y0[no], double ro[no], int maxvertices, int maxfunctions, double *vertices, int *v, FUNCTION *functions, int *f)

Computes the full geometry for a single zenith angle circle. This is the limiting case of AddZenithAngleEllipse for limb-darkened planets or those at full/new phase. Computes the functional form of the circle, its integral, and all points of intersection with the limb of the occultor. Adds the vertices and functions to the running lists.

Parameters
  • zenith_angle: The zenith angle in radians
  • 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
  • maxvertices: Maximum number of vertices in the problem
  • maxfunctions: Maximum number of functions in the problem
  • vertices: The list of vertices
  • v: The number of vertices so far
  • functions: The list of functions
  • f: The number of functions so far

void AddOccultors(double r, int no, double x0[no], double y0[no], double ro[no], int maxvertices, int maxfunctions, double *vertices, int *v, FUNCTION *functions, int *f)

Computes the geometry for all occulting bodies. Adds relevant vertices and functions to the running lists.

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
  • maxvertices: Maximum number of vertices in the problem
  • maxfunctions: Maximum number of functions in the problem
  • vertices: The list of vertices
  • v: The number of vertices so far
  • functions: The list of functions
  • f: The number of functions so far

void AddOcculted(double r, int no, double x0[no], double y0[no], double ro[no], int maxvertices, int maxfunctions, double *vertices, int *v, FUNCTION *functions, int *f)

Computes the geometry for the occulted body. Adds the vertices of intersection with each of the occultors and the functional form of the planet disk to the running lists.

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
  • maxvertices: Maximum number of vertices in the problem
  • maxfunctions: Maximum number of functions in the problem
  • vertices: The list of vertices
  • v: The number of vertices so far
  • functions: The list of functions
  • f: The number of functions so far

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