Procedures

The angle between two vectors (in radians).

Rotate a 3x1 vector in space, given an axis and angle of rotation.

Computes the rotation matrix that corresponds to a rotation about the axis k by an angle theta.

Balance a real general matrix and isolate eigenvalues whenever possible.

Form the eigenvectors of a real general matrix from the eigenvectors of matrix output from BALANC.

Barker time of flight equation

== operator for base_class variables. To be equal, they must be the same type and have the same ID.

/= operator for base_class variables. To be equal, they must be the same type and have the same ID.

Computes the box product (scalar triple product) of the three vectors.

Compute B-plane parameters from position and velocity.

Unit test for bplane_module.

Test of the fmin and zeroin functions.

Convert a c_ptr to a string into a Fortran string.

Calculate current elements z(jd) for planet j from jpl data

Compute B-plane parameters from position and velocity -- alternate version.

Returns the year, month, day, hr, min, sec for the specified Julian date.

Convert Cartesian coordinates to modified equinoctial elements (posigrade formulation).

Function computes the geodetic latitude (phi), longitude (lambda) and height (h) of a point related to an ellipsoid defined by its three semiaxes ax, ay and b (0 < b <= ay <= ax) given Cartesian coordinates Xi, Yi, Zi and tolerance (tol). Latitude and longitude are returned in radians.

Cartesian to geodetic for Triaxial Ellipsoid.

Convert Cartesian (x,y,z) to spherical (r,alpha,beta).

Cartesian to Geodetic I

Cartesian into Geodetic II

Compute the complex quotient of two complex numbers.

Compute the first derivative using a 2-point central difference [-h,h].

Compute the first derivative using a 4-point central difference [-2h,-h,h,2h].

Close the ephemeris.

Close the SPICE ephemeris and unload all the kernels.

Compute the first derivative using the complex-step method. This is Equation 6 from Reference [1].

Unit test for the complex_step module.

Compute $\mu$, the normalized CRTBP parameter. It is equal to $M_2 / (M_1 + M_2)$.

Compute the eigenvalues and, optionally, the eigenvectors of a real general matrix.

Wrapper for Kuga/Carrara method.

Wrapper for Lear method.

Wrapper for Mueller method.

Wrapper for normalized Pines method.

Wrapper for Pines method.

Compute the halo orbit monodromy matrix (which is the state transition matrix propagated for one period) The input should be the result from the halo_to_rv_diffcorr routine.

Compute the CRTBP Jacobi constant, given the state.

Compute the coordinates of the libration points (L1,L2,L3,L4,L5). L1-L3 are computed using Newton's method. L4-L5 are known analytically.

Compute the coordinates of the libration points (L1,L2,L3,L4,L5).

Compute the eigenvalues of the monodromy matrix.

Returns only the real eigenvalues and the associated eigenvectors. Wrapper for compute_eigenvalues_and_eigenvectors.

Compute the new step size using the specific method.

Compute the incoming and/or outgoing v-infinity vectors, given the position and velocity of a hyperbola.

Based on the CONVERT subroutine from [1]. Unnormalizes the C,S coefficients.

Cross product of two 3x1 vectors

Computes the cross product matrix, where: cross(a,b) == matmul(cross_matrix(a),b)

CRTBP derivatives: state only.

CRTBP derivatives: state + state transition matrix.

Unit tests for CRTBP routines.

Cube root of a number (real solution only).

Clohessy-Wiltshire equations for relative motion.

Clohessy-Wiltshire propagation routine.

Derivative of kepde w.r.t de.

Derivative of kepdh w.r.t dh.

Derivative of kepe w.r.t. e

Solution to a*x**3 + 3*b*x - 2c = 0, where a and b**3 + a*c**2 are both non-negative (zero generated, in lieu of infinity, if a = b = 0)

Cube root computed accurately, by incorporating one Newton-Raphson iteration.

Destructor for rk_variable_step_class.

Destructor for rk_class.

This is just a wapper for destroy in geopotential_model.

Destroy a gravity model.

Destructor for stepsize_class.

Solve the "direct" geodetic problem: given the latitude and longitude of one point and the azimuth and distance to a second point, determine the latitude and longitude of that second point. The solution is obtained using the algorithm by Vincenty.

Unit test for the direct and inverse geodetic routines.

Compute the distance between the point X and the line defined by the two points X1 and X2.

Compute the distance between a line segment and a point.

Compute the distance between a point and a polygonal path. Given a point (x0,y0), and a path (x(n),y(n)), the distance to the path is the distance to the closest line segment (x(i),y(i)).

given an m by n matrix a, an n by n nonsingular diagonal matrix d, an m-vector b, and a positive number delta, the problem is to determine the convex combination x of the gauss-newton and scaled gradient directions that minimizes (ax - b) in the least squares sense, subject to the restriction that the euclidean norm of dx be at most delta.

Replacement for the original Minpack routine.

Acceleration due to atmospheric drag.

converts cartesian heliocentric j2000 ecliptic to equatorial

Constructor for a ecliptic_frame

Unit test

Kepler's equation, em = ekepl - (1 - e1)*sin(ekepl), with e1 in range 1 to 0 inclusive, solved accurately (based on ekepl3, but entering e1, not e)

Solve kepler's equation, em = ekepl - e*sin(ekepl), with legendre-based starter and halley iterator (function has also been used under the name eafkep)

Kepler's equation, em = ekepl - e*sin(ekepl) with e in range 0 to 1 inclusive, solved accurately

heliocentric coordinates for orbital plane from elements

Reduce a real general matrix to upper Hessenberg form using stabilized elementary similarity transformations.

Algorithm for two-dimensional conversion from orbital elements to position and velocity.

Algorithm for three-dimensional conversion from orbital elements to position and velocity

Accumulates the stabilized elementary similarity transformations used in the reduction of a real general matrix to upper Hessenberg form by ELMHES.

Similar to emkepl, except input is 1-e.

Accurate computation of ee - e*sin(ee) when (e, ee) is close to (1, 0)

given an n-vector x, this function calculates the euclidean norm of x.

Ephemeris test routine.

converts cartesian heliocentric equatorial to ecliptic

Rotation matrix from J2000 to Mean Ecliptic.

Convert modified equinoctial elements (posigrade formulation) to Cartesian coordinates.

Convert ephemeris time (seconds from J2000 epoch) to Julian date.

Extract the value from the vector and update the index

Extract a vector from the vector and update the index

Convert a Fortran string to a c_ptr to a string. (the C string must already have been allocated to a fixed size)

this subroutine computes a forward-difference approximation to the n by n jacobian matrix associated with a specified problem of n functions in n variables. if the jacobian has a banded form, then function evaluations are saved by only approximating the nonzero terms.

Put the value in the vector and update the index (character version)

Put the vector in the vector and update the index (character version)

Put the value in the vector and update the index

Put the vector in the vector and update the index

The FL factorial function from [1].

An approximation x to the point where f attains a minimum on the interval (ax,bx) is determined.

Compute the first derivative using a forward difference. This is Equation 1 from Reference [1].

Conversion from IJK to LVLH

Compute the transformation matrices to convert IJK to LVLH.

Transform a position (and optionally velocity) vector from IJK to LVLH.

Conversion from IJK to RSW

Compute the transformation matrices to convert IJK to RSW.

Transform a position (and optionally velocity) vector from IJK to RSW.

Conversion from LVLH to IJK

Compute the transformation matrices to convert LVLH to IJK.

Transform a position (and optionally velocity) vector from LVLH to IJK.

Conversion from LVLH to RSW

Transform a position (and optionally velocity) vector from LVLH to RSW.

returns the state of the frame center w.r.t. the frame primary body.

Conversion from RSW to IJK

Compute the transformation matrices to convert RSW to IJK.

Transform a position (and optionally velocity) vector from RSW to IJK.

Conversion from RSW to LVLH

Transform a position (and optionally velocity) vector from RSW to LVLH.

The distance from the center of a celestial body (e.g., the Earth) to a point on the spheroid surface at a specified geodetic latitude.

Geodetic latitude, longitude, and height to Cartesian position vector.

Function computes the Cartesian coordinates given the geodetic latitude (phi), longitude (lambda) and height (h) of a point related to an ellipsoid defined by its three semiaxes ax, ay and b

Geodetic to Cartesian for Triaxial Ellipsoid.

Unit test routine

Compute the gravitational acceleration vector using the Mueller method.

Unit test routine for geopotential_module

This is just a wapper for get_acc in geopotential_model.

rotation matrix for ICRF <-> Mean Ecliptic

rotation matrix for IAU_EARTH <-> ICRF

rotation matrix for IAU_MOON <-> ICRF

rotation matrix for ICRF <-> ICRF

rotation matrix for ROTATING <-> ICRF

rotation matrix for ROTATING_PULSATING <-> ICRF

Obtain the constants from the ephemeris file.

Returns the format statement from a line in a .GEO gravity coefficient file.

Returns a uniform random number x, such that: a <= x < b.

Interface for the ephemeris_module.

Interface for the ephemeris_module.

This subroutine reads the JPL planetary ephemeris and gives the position and velocity of the point ntarg with respect to ncent.

Based on the GRAV subroutine from [1].

Gravitational acceleration due to simplified spherical harmonic expansion (only the J2-J4 terms are used).

Spencer's implementation of the Pines algorithms from [1]

Great circle distance on a spherical body, using the Vincenty algorithm.

Unit test for the halo orbit module.

Compute the state vector from the halo orbit approximation. This will be an approximation of a halo orbit in the CR3BP system, and will need to be corrected to produce a real halo orbit.

Compute the state vector for a halo orbit. This uses the approximation, which is retargeted in the real CR3BP system to produce a periodic orbit.

Heikkinen routine for cartesian to geodetic transformation

For planet np and julian date jd and using using table itbl, return j2000 ecliptic position (au) and velocity (au/yr). in cartesian coordinates (p(1)-p(6)).

computation of an initial step size guess

Horner's method to compute B(x-c) in terms of B(x).

Compute the eigenvalues of a real upper Hessenberg matrix using the QR method.

Compute the eigenvalues and eigenvectors of a real upper Hessenberg matrix using QR method.

Computes a starting step size to be used in solving initial value problems in ordinary differential equations.

The purpose of hybrd is to find a zero of a system of n nonlinear functions in n variables by a modification of the powell hybrid method. the user must provide a subroutine which calculates the functions. the jacobian is then calculated by a forward-difference approximation.

the purpose of hybrd1 is to find a zero of a system of n nonlinear functions in n variables by a modification of the powell hybrid method. this is done by using the more general nonlinear equation solver hybrd. the user must provide a subroutine which calculates the functions. the jacobian is then calculated by a forward-difference approximation.

the purpose of hybrj is to find a zero of a system of n nonlinear functions in n variables by a modification of the powell hybrid method. the user must provide a subroutine which calculates the functions and the jacobian.

the purpose of hybrj1 is to find a zero of a system of n nonlinear functions in n variables by a modification of the powell hybrid method. this is done by using the more general nonlinear equation solver hybrj. the user must provide a subroutine which calculates the functions and the jacobian.

Compute the hyperbolic turning angle from the eccentricity.

Constructor for a iau_earth_rotating_frame

Constructor for a iau_moon_rotating_frame

Returns the rotation matrix for a coordinate transformation from the International Celestial Reference Frame (ICRF) frame to the IAU rotating frame associated with a body. The IAU orientation models use three Euler angles to describe the pole and prime meridian location (ra, dec, and w).

Unit test routine for iau_module.

Constructor for a icrf_frame

Rotation matrix from ICRF to IAU_EARTH.

Rotation matrix from ICRF to IAU_MOON.

Initialize the rk_variable_step_class.

Initialize the rk_class.

Initialize the ephemeris. This routine may be called to load a different ephemeris file. Otherwise, it is called on the first call to get_state, and loads the file specified in the module header.

This is just a wapper for initialize in geopotential_model.

Initialize a SPICE ephemeris by loading the specified kernels.

Main integration routine for the rk_variable_step_class.

Main integration routine for the rk_class.

Event-finding integration routine for the rk_variable_step_class. Integrates until g(t,x)=0, or until t=tf (whichever happens first).

Event-finding integration routine for the rk_class. Integrates until g(t,x)=0, or until t=tf (whichever happens first).

this subroutine differentiates and interpolates a set of chebyshev coefficients to give position and velocity.

INVERSE computes the geodetic azimuth and distance between two points, given their geographic positions.

Convert Julian date to ephemeris time (seconds from J2000 epoch).

Converts Julian date to Modified Julian date.

calendar date to julian date

Returns the Julian date for the specified YEAR, MONTH, DAY, HR, MIN, SEC.

Returns the Julian date for the specified YEAR, MONTH, DAY, HR, MIN, SEC.

Returns the Julian day number (i.e., the Julian date at Greenwich noon) on the specified YEAR, MONTH, and DAY.

Elliptic Kepler's equation written in terms of the eccentric anomaly difference. See Battin, eqn 4.43.

Battin, eqn. 4.64.

Elliptic Kepler's equation

solve kepler's equation ma = ea + ec*sin(ea)

Classical Kepler propagator for elliptical and hyperbolic orbits. Uses Lagrange formulations from Battin & Newton's method.

Kepler propagator based on the Goodyear code with modifications by Stienon and Klumpp.

Kepler propagation using Shepperd's method.

Compute geopotential acceleration using the Kuga/Carrara algorithm. Based on Leg_ForCol_Ac from [1].

Compare the Lambert routines.

given a polygonal line connecting the vertices (x(i),y(i)) (i = 1,...,n) taken in this order. it is assumed that the polygonal path is a loop, where (x(n),y(n)) = (x(1),y(1)) or there is an arc from (x(n),y(n)) to (x(1),y(1)).

Convert the string to lowercase.

Returns a positive number the same magnitude as the input, with only one significant digit.

Compute the cofactors matrix (the transpose of the adjugate matrix).

Matrix determinant of an $n \times n$ matrix (recursive formulation).

Compute the matrix trace (sum of the diagonal elements).

Use maxval(abs(x)) for computing the vector norm.

Rotation matrix from Mean Ecliptic to J2000.

Mean obliquity of the ecliptic, IAU 1980 formula.

Mean obliquity of the ecliptic, IAU 2006 formula.

Converts Modified Julian date to Julian date.

Modified equinoctial elements (posigrade formulation) equations of motion.

Unit tests for the modified_equinoctial_module.

Newton's method for root finding of scalar function f(x)

Use intrinsic norm2(x) for computing the vector norm.

Convert state in km, km/s units to normalized CRTBP state.

Number of (c,s) coefficients for n x m geopotential model Starting with n=2,m=0.

Olson routine for cartesian to geodetic transformation.

heliocentric coordinates j2000 ecliptic plane from orbital plane

Check the orbit for singularities.

Compute the two-body orbital energy.

Compute the two-body orbital period.

Convert orbital elements to position and velocity vectors.

Computes the outer product of the two vectors.

Compute the periapsis and apoapsis position and velocity magnitudes.

Normalized Pines geopotential code.

Print a matrix to the console.

Basic two-body propagator using the Gooding universal element routines.

Algorithm for two-dimensional conversion from position and velocity to orbital elements.

Algorithm for three-dimensional conversion from position and velocity to orbital elements.

this subroutine proceeds from the computed qr factorization of an m by n matrix a to accumulate the m by m orthogonal matrix q from its factored form.

this subroutine uses householder transformations with column pivoting (optional) to compute a qr factorization of the m by n matrix a. that is, qrfac determines an orthogonal matrix q, a permutation matrix p, and an upper trapezoidal matrix r with diagonal elements of nonincreasing magnitude, such that ap = qr. the householder transformation for column k, k = 1,2,...,min(m,n), is of the form

given an m by n matrix a, this subroutine computes aq where q is the product of 2(n - 1) transformations

given an m by n lower trapezoidal matrix s, an m-vector u, and an n-vector v, the problem is to determine an orthogonal matrix q such that

Read the gravity coefficient file. Example file: ftp://ftp.csr.utexas.edu/pub/grav/EGM96.GEO.Z

Unit tests for the relative_motion_module.

Just a test of f_string_to_c_ptr.

Compute the eigenvalues and, optionally, the eigenvectors of a real general matrix.

Take one Runge Kutta 4 integration step: t -> t+h (x -> xf)

Take one Runge Kutta 8 integration step: t -> t+h (x -> xf) This is Formula (8-10) from Reference [1].

Unit test of the rk_module. Integrate a two-body orbit around the Earth.

Unit test of the rk_module. Integrate a two-body orbit around the Earth.

Feagin's RK8(10) method -- a 10th-order method with an embedded 8th-order method.

Returns the order of the rkf108 method.

Feagin's RK12(10) method -- a 12th-order method with an embedded 10th-order method.

Returns the order of the rkf1210 method.

Feagin's RK14(12) - a 14th-order method with an embedded 12th-order method.

Returns the order of the rkf1412 method.

Fehlberg's 7(8) algorithm.

Returns the order of the rkf78 method.

Fehlberg 8(9) method.

Returns the order of the rkf89 method.

Runge Kutta Verner 8(9)

Returns the order of the rkv89 method.

The 3x3 rotation matrix for a rotation about the x, y, or z-axis.

Time derivative of the 3x3 rotation matrix for a rotation about the x, y, or z-axis.

Convert position and velocity vectors to orbital elements.

Returns the state of the from frame center w.r.t. the to frame center, at the specified ephemeris time et.

Set the function to be minimized.

Equation el = shkepl + (g1 - 1)*asinh(shkepl), with g1 in range 0 to 1 inclusive, solved accurately.

Accurate computation of s - (1 - g1)*asinh(s) when (g1, s) is close to (0, 0)

A simple analytical lunar ephemeris model. Returns Lunar cartesian coordinates (mean equator and equinox of epoch J2000).

Solve Lambert's problem using the Arora/Russell method.

Solve Lambert's problem using Gooding's method.

Solve Lambert's problem using Izzo's method.

Numerical solution to polynomial equation using Newton-Raphson method

Special cases for lat/lon/altitude

cartesian to spherical coordinates (angles in radians)

Computes the sphere-of-influence radius of the secondary body.

Computes the sphere-of-influence radius of the secondary body.

Convert spherical (r,alpha,beta) to Cartesian (x,y,z).

Convert the NAIF SPICE ID code to the old one used by the JPL ephemeris. Returns 0 if the body was not found.

Convert the NAIF SPICE ID code to the one used by the standish ephemeris. Returns 0 if the body was not found.

this subroutine breaks a d.p. number into a d.p. integer and a d.p. fractional part.

Test routine for the standish_module routines.

Return the state of the targ body relative to the obs body, in the inertial frame [ICRF].

This subroutine reads and interpolates the JPL planetary ephemeris file.

Unit tests for step size adjustment routines.

Constructor for a stepsize_class.

Stumpff function C(z)

Stumpff function S(z)

Determine which data set to use for highest accuracy for the given julian date.

Third-body (pointmass) gravitational acceleration.

Test routine for the Julian date routines.

Transform a Cartesian state from one reference frame to another at a specified epoch. The from and to reference_frames may each be defined at a different epoch. The et ephemeris time is the time the transformation is to be done, and accounts for the motion of the two frame centers from from%et and to%et to et.

Transformation units test

Constructor for a two_body_rotating_frame

Constructor for a two_body_rotating_pulsating_frame

Unit vector of the cross product of two 3x1 vectors

Time derivative of a unit vector.

Unit vector

Convert normalized CRTBP state to km, km/s units.

Convert the string to uppercase.

The projection of one vector onto another vector.

Project a vector onto a plane.

Unit test routine for the vector_module.

Convert V-infinity magnitude to energy.

Wrap an angle (in rad) from -pi to pi.

Computes the transformation of Cartesian to geodetic coordinates on the surface of the ellipsoid assuming x,y,z are all non-negative Angular coordinates in radians

Determination of the geodetic latitude and longitude

Find a zero of the function $f(x)$ in the given interval $[a_x,b_x]$ to within a tolerance $4 \epsilon |x| + tol$, where $\epsilon$ is the relative machine precision defined as the smallest representable number such that $1.0 + \epsilon > 1.0$.