```````````````````````````````````` ` orblib.num May 21, 2004 ` ```````````````````````````````````` function adb2xyz (adbarv, r, v) ` transform from spherical (adbarv) ` coordinates to cartesian coordinates ` input ` abbarv(1) = geocentric right ascension (radians) ` abbarv(2) = geocentric declination (radians) ` abbarv(3) = conjugate of the flight path angle (radians) ` abbarv(4) = azimuth (radians) ` abbarv(5) = position magnitude ` abbarv(6) = velocity magnitude ` output ` r = eci position vector ` v = eci velocity vector -------------------------------------------------- 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 atmos76 (h, y) ` U.S. Standard 1976 atmosphere model ` linear interpolation - 0 to 1000 kilometers ` input ` h = altitude (kilometers) ` output ` y = data value -------------------------------------------------- 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 coe2eqoe(coe, eqoe) ` convert classical orbital elements to ` equinoctial orbital elements ` input ` coe(1) = semimajor axis (kilometers) ` coe(2) = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` coe(3) = orbital inclination (radians) ` (0 <= inclination <= pi) ` coe(4) = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` coe(5) = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` coe(6) = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) ` output ` eqoe(1) = semimajor axis (kilometers) ` eqoe(2) = e * sin(argument of perigee + raan) ` eqoe(3) = e * cos(argument of perigee + raan) ` eqoe(4) = tan(i/2) * sin(raan) ` eqoe(5) = tan(i/2) * cos(raan) ` eqoe(6) = mean anomaly + raan + argument of perigee -------------------------------------------------- function ecf2eci (gast, recf, vecf, reci, veci) ` ecf-to-eci transformation ` input ` gast = apparent sidereal time (radians) ` (0 <= gast <= 2 pi) ` recf = position vector (kilometers) ` vecf = velocity vector (kilometers/second) ` output ` reci = position vector (kilometers) ` veci = velocity vector (kilometers/second) -------------------------------------------------- function eci2ecf (gast, reci, veci, recf, vecf) ` eci-to-ecf transformation ` input ` gast = apparent sidereal time (radians) ` (0 <= gast <= 2 pi) ` reci = position vector (kilometers) ` veci = velocity vector (kilometers/second) ` output ` recf = position vector (kilometers) ` vecf = velocity vector (kilometers/second) ` global constants ` omega = earth rotation rate (radians/second) -------------------------------------------------- function eci2lvlh (r, v, upeci, uplvlh, pitch, yaw) ` convert eci unit pointing vector to local ` vertical local horizontal coordinate system ` input ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) ` upeci = eci unit pointing vector (non-dimensional) ` output ` uplvlh = lvlh unit pointing vector (non-dimensional) ` pitch = lvlh pitch angle (radians) ` (-pi/2 <= pitch <= +pi/2) ` yaw = lvlh yaw angle (radians) ` (0 <= yaw <= 2 pi) -------------------------------------------------- function eci2eqoe(r, v, eqoe) ` convert eci state vector to ` equinoctial orbital elements ` input ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) ` output ` eqoe[1] = semimajor axis (kilometers) ` eqoe[2] = e * sin(argument of perigee + raan) ` eqoe[3] = e * cos(argument of perigee + raan) ` eqoe[4] = tan(i/2) * sin(raan) ` eqoe[5] = tan(i/2) * cos(raan) ` eqoe[6] = mean anomaly + raan + argument of perigee -------------------------------------------------- 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 eci2rtn1(oev, tm) ` eci-to-rtn transformation matrix ` orbital elements version ` 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) ` output ` tm = eci-to-rtn transformation matrix -------------------------------------------------- function eci2rtn2(r, v, tm) ` eci-to-rtn transformation matrix ` state vector version ` input ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) ` output ` tm = eci-to-rtn transformation matrix -------------------------------------------------- function eqoe2coe (eqoe, coe) ` convert equinoctial orbital elements ` to classical orbital elements ` input ` eqoe(1) = semimajor axis (kilometers) ` eqoe(2) = e * sin(argument of perigee + raan) ` eqoe(3) = e * cos(argument of perigee + raan) ` eqoe(4) = tan(i/2) * sin(raan) ` eqoe(5) = tan(i/2) * cos(raan) ` eqoe(6) = mean anomaly + raan + argument of perigee ` output ` coe(1) = semimajor axis (kilometers) ` coe(2) = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` coe(3) = orbital inclination (radians) ` (0 <= inclination <= pi) ` coe(4) = argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` coe(5) = right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` coe(6) = true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function eqoe2eci(eqoe, r, v) ` convert equinoctial orbital elements ` to eci state vector ` input ` eqoe[1] = semimajor axis (kilometers) ` eqoe[2] = e * sin(argument of perigee + raan) ` eqoe[3] = e * cos(argument of perigee + raan) ` eqoe[4] = tan(i/2) * sin(raan) ` eqoe[5] = tan(i/2) * cos(raan) ` eqoe[6] = mean anomaly + raan + argument of perigee ` output ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) -------------------------------------------------- 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 geodet1 (rmag, dec, alt, lat) ` geodetic latitude and altitude ` series solution ` input ` rmag = geocentric radius (kilometers) ` dec = geocentric declination (radians) ` (+north, -south -pi/2 <= dec <= +pi/2) ` output ` alt = geodetic altitude (kilometers) ` lat = geodetic latitude (radians) ` (+north, -south -pi/2 <= lat <= +pi/2) ` global constants ` req = equatorial radius (kilometers) ` flat = flattening factor (non-dimensional) ` reference: NASA TN D-7522 -------------------------------------------------- function geodet2(rmag, decl, alt, xlat) ` geocentric to geodetic coordinates ` exact solution (Borkowski, 1989) ` input ` rmag = geocentric distance (kilometers) ` decl = geocentric declination (radians) ` output ` alt = geodetic altitude (kilometers) ` xlat = geodetic latitude (radians) -------------------------------------------------- function geodet3 (lat, long, alt, r) ` convert geodetic coordinates to ecf position vector ` input ` lat = geodetic latitude (radians) ` (+north, -south -pi/2 <= lat <= +pi/2) ` long = geodetic longitude (radians) ` alt = geodetic altitude (kilometers) ` output ` r = ecf position vector (kilometers) ` global constants ` req = equatorial radius (kilometers) ` flat = flattening factor (non-dimensional) -------------------------------------------------- function geodet4 (lat, alt, rmag, dec) ` geodetic to geocentric coordinates ` input ` lat = geodetic latitude (radians) ` (+north, -south -pi/2 <= lat <= +pi/2) ` alt = geodetic altitude (kilometers) ` output ` rmag = geocentric position magnitude (kilometers) ` dec = geocentric declination (radians) ` (+north, -south -pi/2 <= dec <= +pi/2) -------------------------------------------------- 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 getsv(r, v) ` request state vector -------------------------------------------------- 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 gravity (t, y, agrav) ` first order equations of orbital motion ` n degree and m order ` gravitational acceleration ` input ` t = simulation time (seconds) ` y = state vector ` output ` agrav = eci gravitational acceleration ` vector (km/sec/sec) -------------------------------------------------- function jacchia (ida, mn, iyr, ihr, xlong, xlat, alt, f10, f10b, gi, dens) ` jacchia 1970 atmosphere model -------------------------------------------------- function j2eqm (t, y, ydot) ` first order equations of orbital motion ` j2 gravitational acceleration ` input ` t = simulation time (seconds) ` y = state vector ` output ` ydot = integration vector -------------------------------------------------- function j4eqm(t, y, ydot) ` first order equations of orbital motion ` includes j2, j3, j4 ` input ` t = simulation time (seconds) ` y = state vector ` output ` ydot = integration vector -------------------------------------------------- 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 kozai1 (iniz, t, oev1, oev2, r, v) ` analytic orbit propagation ` Kozai's method - ECI version ` input ` iniz = initialization flag (1 = initialize, 0 = bypass) ` t = elapsed simulation time (seconds) ` initial orbital elements ` oev1[1] = semimajor axis (kilometers) ` oev1[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev1[3] = orbital inclination (radians) ` (0 <= oev1[3] <= pi) ` oev1(4) = argument of perigee (radians) ` (0 <= oev1(4) <= 2 pi) ` oev1(5) = right ascension of ascending node (radians) ` (0 <= oev1(5) <= 2 pi) ` oev1(6) = mean anomaly (radians) ` (0 <= oev1(6) <= 2 pi) ` output ` final orbital elements and state vector at time = t ` oev2[1] = semimajor axis (kilometers) ` oev2[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev2[3] = orbital inclination (radians) ` (0 <= oev2[3] <= pi) ` oev2(4) = updated argument of perigee (radians) ` (0 <= oev2(4) <= 2 pi) ` oev2(5) = updated right ascension of ascending node (radians) ` (0 <= oev2(5) <= 2 pi) ` oev2(6) = updated mean anomaly (radians) ` (0 <= oev2(6) <= 2 pi) ` oev2(7) = updated true anomaly (radians) ` (0 <= oev2(6) <= 2 pi) ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) -------------------------------------------------- function kozai2 (iniz, t, oev1, oev2, r, v) ` analytic orbit propagation ` Kozai's method - ECF version ` input ` iniz = initialization flag (1 = initialize, 0 = bypass) ` t = elapsed simulation time (seconds) ` initial orbital elements ` oev1[1] = semimajor axis (kilometers) ` oev1[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev1[3] = orbital inclination (radians) ` (0 <= oev1[3] <= pi) ` oev1(4) = argument of perigee (radians) ` (0 <= oev1(4) <= 2 pi) ` oev1(5) = east longitude of ascending node (radians) ` (0 <= oev1(5) <= 2 pi) ` oev1(6) = mean anomaly (radians) ` (0 <= oev1(6) <= 2 pi) ` output ` final orbital elements and state vector at time = t ` oev2[1] = semimajor axis (kilometers) ` oev2[2] = orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oev2[3] = orbital inclination (radians) ` (0 <= oev2[3] <= pi) ` oev2(4) = updated argument of perigee (radians) ` (0 <= oev2(4) <= 2 pi) ` oev2(5) = updated east longitude of ascending node (radians) ` (0 <= oev2(5) <= 2 pi) ` oev2(6) = updated mean anomaly (radians) ` (0 <= oev2(6) <= 2 pi) ` oev2(7) = updated true anomaly (radians) ` (0 <= oev2(6) <= 2 pi) ` r = ecf position vector (kilometers) ` v = ecf velocity vector (kilometers/second) -------------------------------------------------- function kquad(a, b, tol, result, k, icheck) ` one-dimensional quadrature ` method due to T. Patterson with modifications ` by F. Krogh and W. Snyder (ACM Algorithm 699) ` input ` a = lower integration limit ` b = upper integration limit ` tol = convergence tolerance ` output ` result = array of solution(s) ` k = index of result array which best satisfies tol ` icheck = 0 ==> convergence ` = 1 ==> non-convergence -------------------------------------------------- function lambfunc(ri, rf, tof, direct, revmax, statev, nsol) ` solve Lambert's orbital two point boundary value problem ` input ` ri = initial ECI position vector (kilometers) ` rf = final ECI position vector (kilometers) ` tof = time of flight (seconds) ` direct = transfer direction (1 = posigrade, -1 = retrograde) ` revmax = maximum number of complete orbits ` output ` nsol = number of solutions ` statev = matrix of solutions -------------------------------------------------- function interp1(x, y, xval, fval) ` linear interpolation subroutine ` y = f(x) ` input ` x = vector of x data (n rows) ` y = vector of y data (n rows) ` xval = x argument ` output ` fval = interpolated function value at xval -------------------------------------------------- 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 mean2osc (oemean, oeosc) ` convert mean classical orbital elements to ` osculating classical orbital elements ` input ` oemean[1] = mean semimajor axis (kilometers) ` oemean[2] = mean orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oemean[3] = mean orbital inclination (radians) ` (0 <= inclination <= pi) ` oemean[4] = mean argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oemean[5] = mean right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oemean[6] = mean true anomaly (radians) ` (0 <= true anomaly <= 2 pi) ` output ` oeosc[1] = osculating semimajor axis (kilometers) ` oeosc[2] = osculating orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oeosc[3] = osculating orbital inclination (radians) ` (0 <= inclination <= pi) ` oeosc[4] = osculating argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oeosc[5] = osculating right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oeosc[6] = osculating true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- 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 nym4(n, tp, dt, rs, vs, rf, vf) ` solve a system of second order differential equations ` fourth order Nystrom method (fixed step size) ` input ` n = number of equations in user-defined ` system of differential equations ` tp = current simulation time ` dt = integration step size ` rs = position vector at initial time ` vs = velocity vector at initial time ` output ` rf = position vector at time = tp + dt ` vf = velocity vector at time = tp + dt -------------------------------------------------- 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 orbital 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 orbital 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 osc2mean (oeosc, oemean) ` convert osculating classical orbital elements ` to mean classical orbital elements ` input ` oeosc[1] = osculating semimajor axis (kilometers) ` oeosc[2] = osculating orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oeosc[3] = osculating orbital inclination (radians) ` (0 <= inclination <= pi) ` oeosc[4] = osculating argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oeosc[5] = osculating right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oeosc[6] = osculating true anomaly (radians) ` (0 <= true anomaly <= 2 pi) ` output ` oemean[1] = mean semimajor axis (kilometers) ` oemean[2] = mean orbital eccentricity (non-dimensional) ` (0 <= eccentricity < 1) ` oemean[3] = mean orbital inclination (radians) ` (0 <= inclination <= pi) ` oemean[4] = mean argument of perigee (radians) ` (0 <= argument of perigee <= 2 pi) ` oemean[5] = mean right ascension of ascending node (radians) ` (0 <= raan <= 2 pi) ` oemean[6] = mean true anomaly (radians) ` (0 <= true anomaly <= 2 pi) -------------------------------------------------- function planet(ip, jdate, r, v) ` planetary ephemeris ` input ` ip = planet index (1 = mercury, 2 = venus, etc.) ` jdate = julian day ` output ` r = heliocentric position vector (kilometers) ` v = heliocentric velocity vector (kilometers/second) ` Note: coordinates are with respect to the mean ` ecliptic and equinox of date. -------------------------------------------------- function polyroot(coef, npoly, rootr, rooti, nflg) ` real and complex roots of a polynomial ` f(x) = coef(npoly+1) * x ^ npoly + coef(npoly) * x ^ (npoly-1) ` + coef(npoly-1) * x ^ (npoly-2) + ... + coef(1) ` input ` coef = array of polynomial coefficents ` npoly = order of the polynomial ` (1 <= order <= 36) ` output ` rootr = array of real parts of roots ` rooti = array of imaginary parts of roots ` nflg = error flag ` (nflg = 0 ==> no error) ` (nflg = 1 ==> iterations > 500) -------------------------------------------------- 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 svprint(r, v) ` print state vector ` r in kilometers ` v in kilometers/second -------------------------------------------------- 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 read76(ad76) ` read U.S. Standard 76 atmosphere data file ` output ` ad76 = atmospheric density data array -------------------------------------------------- function readgrav(fname, ccoef, scoef) ` read gravity model data file ` input ` fname = name of data file ` output ` ccoef, scoef = gravity model data -------------------------------------------------- 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 readtle(satstr, tlename) ` read NORAD two line element file ` and extract orbital information ` input ` satstr = name of satellite ` tlename = name of TLE data file ` output ` fid = file id ` NOTE: all TLE output via global -------------------------------------------------- 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 sgp4 (tsince, r, v) ` SGP4 orbit propagation ` input ` tsince = time since initialization (minutes) ` output ` r = position vector (kilometers) ` v = velocity vector (km/sec) -------------------------------------------------- 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 stm2 (mu, tau, ri, vi, rf, vf, stm) ` two body state transition matrix ` Shepperd's method ` input ` mu = gravitational constant ` tau = propagation time interval (seconds) ` ri = initial eci position vector (kilometers) ` vi = initial eci velocity vector (km/sec) ` output ` rf = final eci position vector (kilometers) ` vf = final eci velocity vector (km/sec) ` stm = state transition matrix -------------------------------------------------- function stumpff(e2, c1, c2, c3) ` Stumpff functions for argument E^2 ` input ` e2 = eccentric anomaly squared ` output ` c1, c2, c3 = Stumpff functions -------------------------------------------------- function tunion(a, na, b, nb, c, nc) ` determines the logical union of ` two sets of time intervals ` input ` a = array of endpoints of intervals in set a ` na = number of intervals in set a ` b = array of endpoints of intervals in set b ` nb = number of intervals in set b ` output ` c = array of endpoints of intervals in union set c ` nc = number of intervals in union set c -------------------------------------------------- function twobody2 (mu, tau, ri, vi, rf, vf) ` solution of the two body initial value problem ` Shepperd's method ` input ` tau = propagation time interval (seconds) ` ri = initial eci position vector (kilometers) ` vi = initial eci velocity vector (km/sec) ` output ` rf = final eci position vector (kilometers) ` vf = final eci velocity vector (km/sec) -------------------------------------------------- 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 utvector (r, v, pitch, yaw, utv) ` eci unit pointing vector ` input ` r = eci position vector (kilometers) ` v = eci velocity vector (kilometers/second) ` pitch = rtn pitch angle (radians) ` yaw = rtn yaw angle (radians) ` output ` utv = eci unit pointing 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 -------------------------------------------------- function xyz2adb (r, v, adbarv) ` transform from cartesian coordinates ` to spherical (adbarv) coordinates ` input ` r = eci position vector ` v = eci velocity vector ` output ` abbarv(1) = geocentric right ascension (radians) ` abbarv(2) = geocentric declination (radians) ` abbarv(3) = conjugate of the flight path angle (radians) ` abbarv(4) = azimuth (radians) ` abbarv(5) = position magnitude ` abbarv(6) = velocity magnitude --------------------------------------------------