function atan3 (a, b) ` four quadrant inverse tangent ` input ` a = sine of angle ` b = cosine of angle ` output ` angle (radians; 0 =< c <= 2 * pi) -------------------------------------------------- function brent (x1, x2, rtol, xroot, froot) ` solve for a single real root of a nonlinear equation ` Brent's method ` input ` x1 = lower bound of search interval ` x2 = upper bound of search interval ` rtol = algorithm convergence criterion ` output ` xroot = real root of f(x) = 0 ` froot = function value at f(x) = 0 -------------------------------------------------- function broot (x1in, x2in, factor, dxmax, x1out, x2out) ` bracket a single root of a nonlinear equation ` input ` x1in = initial guess for first bracketing x value ` x2in = initial guess for second bracketing x value ` factor = acceleration factor (non-dimensional) ` dxmax = rectification interval ` output ` x1out = final value for first bracketing x value ` x2out = final value for second bracketing x value -------------------------------------------------- function eci2orb (mu, r, v, oev) ` convert eci state vector to six ` classical orbital elements via ` equinoctial elements ` input ` mu = gravitational constant (km^3/sec^2) ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) ` output ` oev[1] = semimajor axis (kilometers) ` oev[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev[3] = orbital inclination (radians) ` (0 <= inclination <= pi) ` oev[4] = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oev[5] = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oev[6] = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function events (ti, tf, topt) ` display circumstances of celestial events ` input ` ti = initial simulation time ` tf = final simulation time ` topt = extrema time -------------------------------------------------- function gdate (jdate, month, day, year) ` convert Julian date to Gregorian (calendar) date ` input ` jdate = julian day ` output ` month = calendar month [1 - 12] ` day = calendar day [1 - 31] ` year = calendar year [yyyy] ` note: day may include fractional part -------------------------------------------------- function getdate(cdate, month, day, year) ` interactive calendar date request ` input ` cdate = calendar date string ` output ` month = calendar month ` day = calendar day ` year = calendar year -------------------------------------------------- function gettime(utstr, hours, minutes, seconds) ` interactive time request ` input ` utstr = universal time string ` output ` hours = time in hours ` minutes = time in minutes ` seconds = time in seconds -------------------------------------------------- function getet(etstr, hours, minutes, seconds) ` interactive ephemeris time request ` input ` etstr = ephemeris time string ` output ` hours = time in hours ` minutes = time in minutes ` seconds = time in seconds -------------------------------------------------- function getobs(obslatstr, obslonstr, obslat, obslong, obsalt) ` interactive request of observer coordinates ` input ` obslatstr = observer latitude string ` obslonstr = observer longitude string ` output ` obslat = geographic latitude (radians) ` obslong = geographic longitude (radians) ` obsalt = geodetic altitude (kilometers) -------------------------------------------------- function getoe(ioev, oev) ` interactive request of classical orbital elements ` NOTE: all angular elements are returned in radians ` input ` ioev = request array (1 = yes, 0 = no) ` output ` oev[1] = semimajor ` oev[2] = orbital eccentricity ` oev[3] = orbital inclination ` oev[4] = argument of perigee ` oev[5] = right ascension of the ascending node ` oev[6] = true anomaly -------------------------------------------------- function julian (month, day, year, jdate) ` Julian date ` Input ` month = calendar month [1 - 12] ` day = calendar day [1 - 31] ` year = calendar year [yyyy] ` Output ` jdate = Julian date ` special notes ` (1) calendar year must include all digits ` (2) will report October 5, 1582 to October 14, 1582 ` as invalid calendar dates and stop -------------------------------------------------- function kepler1 (manom, ecc, eanom, tanom) ` solve Kepler's equation for circular, ` elliptic and hyperbolic orbits ` Danby's method ` input ` manom = mean anomaly (radians) ` ecc = orbital eccentricity (non-dimensional) ` output ` eanom = eccentric anomaly (radians) ` tanom = true anomaly (radians) -------------------------------------------------- function kepler2(t, q, ecc, r, tanom) ` Kepler's equation for heliocentric ` parabolic and near-parabolic orbits ` input ` t = time relative to perihelion passage (days) ` q = perihelion radius (AU) ` ecc = orbital eccentricity (non-dimensional) ` output ` r = heliocentric distance (AU) ` tanom = true anomaly (radians) -------------------------------------------------- function mdofunc (n, x, f) ` mdopt user-defined function ` input ` n = number of variables ` x = function argument vector ` output ` f = scalar value of objective function at x -------------------------------------------------- function mdopt(method, gtype, n, eps, maxiter, iflag, niter, f, x) ` multivariable minimization ` numerical gradient method ` input ` method = method of solution ` 1 = conjugate gradient ` 2 = quasi-newton ` gtype = type of gradient calculation ` 1 = Ridder's method ` 2 = central differences ` n = number of variables ` eps = convergence criterion ` maxiter = maximum number of iterations ` x = initial guess for solution vector ` output ` niter = number of algorithm iterations ` f = final objective function value ` x = final solution vector ` iflag = diagnostic flag ` 0 = converged ` 1 = maximum number of function evaluations ` 2 = linear search failure ` 3 = search vector failure -------------------------------------------------- function minima (a, b, tolm, xmin, fmin) ` one-dimensional minimization function ` Brent's method ` input ` a = initial x search value ` b = final x search value ` tolm = convergence criterion ` output ` xmin = minimum x value ` fmin = minimum function value -------------------------------------------------- function modulo (x) ` modulo 2 pi function ` input ` x = argument (radians) ` output ` x modulo 2 pi -------------------------------------------------- function numgrad1 (n, x, fx, gx) ` numerical gradient - Ridder's method ` input ` n = number of variables ` x = vector of independent variables ` output ` fx = scalar function value at x ` gx = gradient values at x -------------------------------------------------- function numgrad2 (n, x, fx, gx) ` numerical gradient - central differences ` input ` n = number of variables ` x = vector of independent variables ` output ` fx = scalar function value at x ` gx = gradient values at x -------------------------------------------------- function objfunc(x, fx) ` "dummy" objective function ` input ` x = function argument ` output ` fx = function value at x -------------------------------------------------- function observer(jed, rsite) ` observer eci position vector ` input ` jed = julian ephemeris date ` output ` rsite = eci position vector -------------------------------------------------- function oeprint(mu, oev) ` print six classical orbital elements ` (orbital period in minutes) ` input ` oev[1] = semimajor axis (kilometers) ` oev[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev[3] = orbital inclination (radians) ` (0 <= inclination <= pi) ` oev[4] = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oev[5] = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oev[6] = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function oeprint1(mu, oev) ` print six classical orbital elements ` (orbital period in days) ` input ` oev[1] = semimajor axis (kilometers) ` oev[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev[3] = orbital inclination (radians) ` (0 <= inclination <= pi) ` oev[4] = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oev[5] = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oev[6] = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function oeprint2(mu, oev) ` print six classical orbital elements ` (orbital period in years, sma in au) ` input ` oev[1] = semimajor axis (kilometers) ` oev[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev[3] = orbital inclination (radians) ` (0 <= inclination <= pi) ` oev[4] = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oev[5] = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oev[6] = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function oevent1 (ti, tf, dt, dtsml) ` predict "extrema" type celestial events ` input ` ti = initial simulation time ` tf = final simulation time ` dt = step size used for bounding minima ` dtsml = small step size used to determine whether ` the function is increasing or decreasing ` output ` via function events -------------------------------------------------- function oevent2 (ti, tf, dt, dtsml) ` predict root-finding celestial events ` input ` ti = initial simulation time ` tf = final simulation time ` dt = step size used for bounding minima ` dtsml = small step size used to determine whether ` the function is increasing or decreasing -------------------------------------------------- function orb2eci(mu, oev, r, v) ` convert classical orbital elements to eci state vector ` input ` mu = gravitational constant (km^3/sec^2) ` oev[1] = semimajor axis (kilometers) ` oev[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev[3] = orbital inclination (radians) ` (0 <= inclination <= pi) ` oev[4] = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oev[5] = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oev[6] = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) ` output ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) -------------------------------------------------- function orbeqm (time, y, ydot) ` first order equations of orbital motion ` cowell's method ` input ` time = simulation time (seconds) ` y = state vector ` output ` ydot = integration vector -------------------------------------------------- function prtdate(month, day, year) ` print calendar date ` input ` imonth = calendar month ` iday = calendar day ` iyear = calendar year -------------------------------------------------- function prthms(angarg) ` print angular element in h m s format ` input ` angarg = anglular element (hours) -------------------------------------------------- function prttime(thours) ` print time in h m s format ` input ` thours = time (hours) -------------------------------------------------- function prtrasc(rasc) ` print right ascension in h m s format ` input ` rasc = right ascension (hours) -------------------------------------------------- function prtangle(angarg) ` print angle in dms format ` input ` angarg = angular argument (degrees) -------------------------------------------------- function r2r (x) ` revolutions to radians function ` input ` x = argument (revolutions; 0 <= x <= 1) ` output ` y = equivalent x (radians; 0 <= y <= 2 pi) -------------------------------------------------- function readoe2(fname, oev, jdpp) ` read heliocentric orbital elements file ` input ` fname = name of orbital elements data file ` output ` oev[1] = perihelion radius (AU) ` oev[2] = orbital eccentricity (non-dimensional) ` oev[3] = orbital inclination (radians) ` oev[4] = argument of perihelion (radians) ` oev[5] = right ascension of ascending node (radians) ` jdpp = julian date of perihelion passage -------------------------------------------------- function rkf78 (neq, ti, tf, h, tetol, x, xout) ` solve first order system of differential equations ` Runge-Kutta-Fehlberg 7(8) method ` input ` neq = number of differential equations ` ti = initial simulation time ` tf = final simulation time ` h = initial guess for integration step size ` tetol = truncation error tolerance (non-dimensional) ` x = integration vector at time = ti ` output ` xout = integration vector at time = tf -------------------------------------------------- function spcoef (x, y, n, c1, c2, c3) ` calculate spline coefficients ` input ` x = x data array ` y = y data array ` n = number of data points in arrays ` output ` c1, c2, c3 = array of spline coefficients -------------------------------------------------- function spdata (x, y, ndata, nsp, xk, yk, npts) ` sample and load data arrays for spline operations ` input ` x = x data array ` y = y data array ` ndata = number of data points in data arrays ` nsp = data points sampling interval ` output ` xk = sampled x data array ` yk = sampled y data array ` npts = number of data points in sampled arrays -------------------------------------------------- function spval (u, x, y, n, c1, c2, c3, sx) ` evaluate spline fit ` input ` u = x value to be fitted ` x = x data array ` y = y data array ` n = number of data points in arrays ` c1, c2, c3 = spline fit coefficient arrays ` output ` sx = spline fit at u value -------------------------------------------------- function snle(n, x, fvec, tol, info) ` solution of a systems of nonlinear equations ` input ` n = number of equations in nonlinear system ` x = initial guess for solution vector ` tol = solution tolerance ` output ` x = final solution vector ` fvec = nonlinear system evaluated at x ` info = solution information flag ` = 1 ==> improper input parameters ` = 2 ==> number of calls to snlefunc has reached ` or exceeded 200 * (n + 1) ` = 3 ==> tol is too small. no further improvement ` in the approximate solution x is possible ` = 4 ==> iteration is not making good progress -------------------------------------------------- function snlefunc(n, x, fvec, iflag) ` user-defined system of nonlinear equations ` input ` n = number of equations in system ` x = current argument vector ` output ` fvec = system of nonlinear equations evaluated at x ` iflag = indicator flag ` < 0 ==> user termination ` >= 0 ==> no action -------------------------------------------------- function stumpff(e2, c1, c2, c3) ` Stumpff functions for argument E^2 ` input ` e2 = eccentric anomaly squared ` output ` c1, c2, c3 = Stumpff functions -------------------------------------------------- function twobody4 (mu, tau, svi, svf) ` solution of the two body initial value problem ` Shepperd's method ` saves the value of u between calls ` input ` mu = gravitational constant ` tau = propagation time interval (seconds) ` svi = initial eci state vector ` output ` svf = final eci state vector -------------------------------------------------- function uvector(x, vhat) ` unit vector ` input ` x = 3-component vector ` output ` vhat = 3-component unit vector -------------------------------------------------- function vcross(a, b, c) ` vector cross product ` input ` a = 3-component vector ` b = 3-component vector ` output ` c = 3-component vector: a cross b -------------------------------------------------- function vdot(x1, x2) ` vector dot product ` input ` x1 = 3-component vector ` x2 = 3-component vector ` output ` dot = vector dot product x1 dot x2 --------------------------------------------------