Home

The User Manual can be found here. - CCAR

1. Orbit integrators or propagators are a way of solving the ordinary differential equation ODE X t F X t t 1 where X is a state vector t is time and X t is the derivative of X with respect to time With an initial state X to and a method of computing F X to to future states X t can be obtained The method of computing the future states is called an initial value solver for the ODE There are eight ODE solvers in TurboProp They are rk4fix rk4fix_dstm rk 45 rk78_dstm rk78 rk78_dstm rk4sym and rk6sym The interface with these solvers varies with the calling software to exploit advantages of the given language A description of these integrators can be called from MATLAB by typing help solvername where solvername is the name of the Turbo Prop integrator Python help information can be accessed by importing the integrator module from the Python interpreter and using the help function The _dstm solvers are useful for sequential processors in orbit determination For example when using a Kalman filter it is much faster to perform one function call to integrate a whole trajectory with automatic initialization of the state transition matrix STM The slower alternative would be to call the integrator in chunks manually reinitializing the STM after each time Although the integrators are very fast there is a performance penalty associated with the function call The interface between the scripting language and the C code r
2. uk f len ODE Function Descriptions lt lt res o 3 9 3 99 o ORC A 39541 3 2 2 2 4 3 5 2 4 24 3 5 2 6 24 7 3 2 5 3 28 32 10 324511 Tal A E 5 2 14 344 12 1 3 2 10 S DT 5 2 18 5 2 19 5 2 20 3 41 DL 3 2 23 5 2 24 ER 2 4 25 Asc IOPU dG 2d RRR cx eom We OR XR SP RRS HERP SRG SESE SEES IWOGBOOUOSER zd xh ox eggs ERE xs Ge ow x eom AUN UR EE eg EPO s Sh cnaRoe rc ve Re BOG edd TROY grow BEB 25v Roe ee wp SER EeeWueus TWOBOdY ESO DRM r enra n Row woe ouo GA deum de Rede des TwoBody Feo Grav SEM sy eid o See eda Gee ed add tas WB GTAV URE 205202 A Se Be eB Be eR a Be E cR d IWOBOOF ALA Grav uk sasear ansa amak na aaea INTE 2x0 nue 65 Gee a a a a Se id LIRTBP EE Lon 20609 2x08 mox Gem oum o me amp Vnus S vie Se e aa ie a CEPS D A E Ede ARD ee A Rte roms ur a er POR TEP SCAD o Loose e wok are eO RO UE mU EOR Em GE oe DE A INE GEM 0906909 5 ew de xeu E ee Note deque RU de DE GIDESVOSLENL ceca teen Shaheen ee es eb wk X cR OR COR S DEORE GTST x0 ok edd E eae Aeg REC RI RO ae DERE eV SEM e eke EDR AA SEE A wed DEPP a ey ea ek RE ERE ER BRE xe qe xeu Se qe RR XU WU Se DE Aire og Pers EM 2 3 6 eG GGA eS om f vom few ow EE EE a DE rag QRP OHOWV reyre GACY GO SOE WE SoG Be Sey DE drag SRP rey St rra wimow w meo ae ee RR we PO XLIEMGSSES 221423949 x Y xo ed x ek wv yx wk rob x xg 1 TURBOPROP 6 Auxiliary Functions 6
3. w We 3 2 74 w3 y 3 3 We is the rotation rate of the Earth and it is multiplied by the third row of the y matrix used to convert a position from inertial to fixed see Section 5 1 6 Applying the cross product gives Vas T Ww9z w3y Vay Y wt Wz 75 Vaz 2 WY wT For DE_drag type functions we is computed automatically and is not supplied by the user 5 1 9 State Transition Matrix stm In stm type derivative functions the state vector is augmented with a state transition matrix STM When the state vector includes a state transition matrix STM it should be added on to the end of the state vector in column order As an example assignment in MATLAB with a six element state and a STM you could type yO x y z dx dy dz this assigns the state y0 6 1 6 36 reshape eye 6 36 1 this adds on the STM State deviations at time t may be mapped back to the initial epoch t with the state transition matrix A sample state transition matrix for two dimensions of position and velocity would be Ox Ox Or Ox Oy Oy Oy Oy Ox Of Oc Of S Ox OYR Oxy OYx Oy Oy Og OY OL OYR Otk OYx P t tk 37 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS The state transition matrix is integrated along with the state using the following relation _ OF X 0 P t tr A t t ti where A t OX gt and is of the form 78 Ox Oy Oz OZ Ox Oy Oz AS EC THE
4. 5 y velocity km s 6 Z velocity km s 7 Cp drag coefficient y0 8 cp reflectance y0 9 72 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields 0 2 yO yO yO yO yO 0 0 0 0 0 extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras degord Section 4 15 6 extras degordstm Section 4 15 7 extras gravfile Section 4 15 8 extras Am m kg Section 4 15 9 extras reftime JED Section 4 15 10 extras atmos kg m km km Section 4 15 12 extras mu km s optional Section 4 15 1 extras radius km optional Section 4 15 3 extras EOP 15 vector Optional when extras center 11 Section 4 15 11 5 2 27 PointMasses PointMasses is a model of the gravity field using the point mass model The state y0 isa 6 x 1 vector u yO 1 x position y0 2 y position y0 3 z position y0 4 x velocity y0 5 y velocity y0 6 z velocity The origin of the coordinate frame is at the center of mass of the point masses The ext ras structure should have all of the following fields
5. The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras Am m kg Section 4 15 9 extras reftime JED Section 4 15 10 extras mu km s optional Section 4 15 1 5 2 19 DE grav DE grav is a model of the solar system using the JPL DE ephemerides and a spherical harmonic gravity model for one of the planets The state y0 is a 6 x 1 vector yO 1 az position km 2 y position km 3 z position km 4 qa velocity km s 5 y velocity km s y0 6 z velocity km s The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields yO yO yO yO extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras reftime JED Section 4 15 10 extras mu km s optional Section 4 15 1 extras radius km optional Section 4 15 3 extras EOP 15 vector Optional when extras center 11 Section 4 15 11 54 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 2 00 DE grav stm DE_grav_stm is a model of the solar
6. TurboProp Image credit NASA Glenn Research Center TurboProp Version 4 0 21 May 2009 Keric Hill and Brandon Jones Project Geryon Colorado Center for Astrodynamics Research University of Colorado Boulder Geryon a warrior in Greek mythology with three bodies was only defeated when Hercules shot him with a poison arrow The aim of the Project Geryon is to conquer the Three body Problem and use its unique astrodynamics to create new possibilities for the exploration of space http ccar colorado edu geryon CONTENTS CONTENTS Contents 1 TurboProp 4 2 Installation 6 2 1 MATLAB Installation J 22 Python Installation III 8 3 Python Interface Definition 9 4 Integrators 10 Wo a EE EEE HS HOw E E E 6 Oe ESSE ESSE 11 de ROTI Oe Lih no See Behe Sed a a da ds 11 Won EEI e E 11 4X ROR Se circadiano 11 do BRR ok tee eee Siw Soke Dae Sots See Sols bode a He bile xo MSS ESE 11 Lo ELA eek Ae Ge ox cv xv AA We BAO WWE eus 12 LT EISI TTL 12 Lo GE SUL a e ii E AA em A quos 12 4 9 MATLAB Integrator Call Syntax 2c 52s RRR ER EER ER ESSE EE 12 4 10 Python Integrator Call Syntax 2 2 ke eee RR 13 LARES lus do o 2 a ee A A AR RAR A AAA RS AA 14 4 12 tspan AIPM csi a ama AAA ee Re a 14 E Ki A ee s eres ie ap OR ae He ROR Re Bee Fee e al a 15 14 CORI ONS TI es es b sos ao e OH RO AAA A A Rae x 15 tD ercran api as ao s eo ee ee BO A Re PS 15 NS See RAS I eran Rw 3 4 4 9 44 c9 3 3 3e Ae Sod Xd eg ES SS 15 15 2 extras
7. mexmac Windows mexw32 or d11 The MEX functions should be reinstalled on each machine on which you wish to use them even if they use the same operating system They are not usually compatible when copied and pasted from one machine to another In Windows the installation script will automatically add the proper directories to the end of the MATLAB path for you so you can call the TurboProp functions just like MATLAB s own built in functions UNIX Macintosh or Linux users will have to save the TurboProp directories in the path definition file or startup file themselves to make it permanent Make sure that the He1p directory is last Another option is to run the TurboPath m script at every startup to add the TurboProp paths to MATLAB for that session All ephemeris files are visible in the MATLAB path so to retrieve the file name and path from any directory use the function which in MATLAB For example to find the DEA05 file you can type extras ephemfile which DE 405 extras ephenmfile will now contain a string with the entire path to the DE405 file 2 2 Python Installation The Python installation is similar to the MATLAB directions provided previously Run the InstallScript py from the Python or shell command line to compile the necessary C code into Python libraries Note there are some dependencies which are summarized below e Simplified Wrapper and Interface Generator SWIG http www swig
8. 79 Ox Oy Oz d M a Ox Oz Although it isn t shown the parameters Cp and cg would also be included in the state transition matrix if they were in the state vector The following sections describe how these partials are computed Central Body STM To include the two body acceleration in the variational equations used to create matrix A the partials of the two body portion of r with respect to r are computed in this manner OF 2b0dy _ He _ H O ody _ TV OF2body _ 02 Ox r5 r3 Oy r5 Oz r5 OY2body ghty OYobody QI H OYavody 3 Hz 80 Ox r5 Oy r5 r3 Oz r5 i O obody phi OZ2body 2 3 FYZ O obody E uz m H Ox r5 Oy r5 Oz r5 r3 where x y and z are the inertial position vector from the center of mass of the body to the spacecraft and r is the magnitude of that vector y is the gravitational parameter of the central body 7 Y and 2 are the accelerations in the x y and z directions respectively Third Body STM To include the third body acceleration in the variational equations used to create matrix A the partials of the third body portion of r with respect to r are computed in this manner OF 3body HaX u3 OX 3boay H3X3Y3 OX3body 3 23 23 ae ee 3 ue uro cm 0x3 r3 r3 Oys r3 0z3 r3 O sbay a H3T3Y3 Oisvody a H3Y3 H3 O sbody ag HB 323 BA 3 3 5 3 81 0x3 T3 Oys T3 r3 023 r3 O 3body H3323 O 3body H3Y323 OZsboay H325 H3 vm ons
9. Input y0 is a vector either a row or a column with the state parameters In Python yO is an array of floats y0 should describe the state at the time given in tspan 1 When the state vector includes a state transition matrix STM it should be added on to the end of the state vector in column order As an example assignment in MATLAB with a six element state and a STM you could type yO x y z dx dy dz this assigns the state yO 6 1 6 36 reshape eye 6 36 1 this adds on the STM In Python the assignment is similar but uses the syntax for assigning variables to an array How ever for definition of the STM you could use for i in xrange 0 6 for j in xrange 0 6 if i j y0 append 1 0 else yO append 0 0 4 14 options Input options is in MATLAB a vector either a row or a column of numbers that are used by the inte grator In Python it is an array of floats options can include integration tolerances step sizes and other parameters used by the integration function See sections 4 1 through 4 6 for details on what to include in the options vector 4 15 extras Input In MATLAB extras is a structure data type passed on to the ODE function or derivative function This variable is used to pass on constants that are needed to compute the derivative of the state but that aren t estimated The fields in extras are described in the following sections If a field is included in extras that is not needed by
10. J B and B E Schutz Recursion Formulas of Legendre Functions for Use with Nonsingu lar Geopotential Models Journal of Guidance Volume 11 1 Jan Feb 1988 pp 31 38 Mohr P J B N Taylor and D B Newell CODATA Recommended Values of the Fundamental Phys ical Constants 2006 Technical report National Institute of Standards and Technology Gaithers burg Maryland 20899 USA December 2007 Murray C D and S F Dermott Solar System Dynamics Cambridge University Press Cambridge UK 1999 Newhall X X and J G Williams Estimation of the Lunar Physical Librations Celestial Mechanics and Dynamical Astronomy Volume 66 1997 pp 21 30 NGA Department of Defense World Geodetic System 1984 Its Definition and Relationships With Local Geodetic Systems Technical Report TR 8350 2 National Geospatial Intelligence Agency January 2000 Standish E M JPL Planetary and Lunar Ephemerides DE405 LE405 Jet Propulsion Laboratory Interoffice Memorandum IOM 312F 98 048 Aug 26 1998 Standish E M X X Newhall J G Williams and W M Folkner JPL Planetary and Lunar Ephemerides DE403 LE403 Jet Propulsion Laboratory Interoffice Memorandum IOM 314 10 127 May 22 1995 Tapley B J Ries S Bettadpur D Chambers M Cheng F Condi B Gunter Z Kang P Nagel R Pastor T Pekker S Poole and F Wang GGM02 An Improved Earth Gravity Field Model from GRACE Journal
11. 1 vector yO 1 ac position y0 2 y position y0 3 z position y0 4 Za velocity y0 5 y velocity y0 6 z velocity The origin of the coordinate frame is at the center of mass of the primary body extras mu must be included and should contain u the gravitational parameter of the primary body It is up to the user to make sure the units of y0 t span and the gravitational parameter are consistent 5 2 2 TwoBody stm This ODE function is a model of the two body problem where a spacecraft orbits a single massive central body The state y0 is a 42 x 1 vector The first 6 values are yO 1 m position y0 2 y position y0 3 z position y0 4 Za velocity y0 5 y velocity y0 6 z velocity y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is at the center of mass of the primary body extras mu must be included and should contain p the gravitational parameter of the primary body It is up to the user to make sure the units of y0 t span and the gravitational parameter are consistent 45 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 2 3 TwoBody grav TwoBody grav is a model of the two body problem where a spacecraft orbits a single massive central body The central body also has a spherical harmonic gravity field The state yO isa6 x 1 vector yO 1 m position y0 2 y position y0 3 z position y0 4 az velocity y0 5 y velocity y0 6
12. 26712767857796x108 1 26712767857796 x 108 Saturn 3 79406260611373x10 3 79406260611373 x10 Uranus 5 79454900707188 x 10 5 79454900707188 x 109 Neptune 6 83653406387926x10 6 83653406387926 x 10 Pluto 9 81600887707005 x10 9 81600887707005 x10 Sun 1 32712440017987 x10 1 32712440017987 x 10H Earth 3 98600435608103 x10 3 98600432896939 x10 Moon 4 90279910787977x10 4 90280058214776 x10 Units are km s 5 1 4 Mass Concentration Gravity Field PointMasses The mass concentration model is a gravity field where the total acceleration is determined by a collec tion of point masses each with a given mass and location The net acceleration is miri Tpointmass G D nm 20 where m is the mass of th point and r is the vector from the point mass to the satellite r is the magnitude of the r vector Note this is essentially the sum of the two body accelerations due to each point mass In TurboProp G 6 67428 x 107 km kg sec 21 This value was selected based on Mohr et al 2007 Note based on the units of this value all input values must be in km and seconds The definition of the point masses sometimes referred to as mascons is provided to the integrator in the extras gravfile option The file format is rather straightforward with one header line and a collection of data lines The format of the header and data are provided in Tables 2 and 3 respectively The given mass for each poin
13. 59 6 AUXILIARY FUNCTIONS extras gravfile Section 4 15 8 extras EOP 2 vector Section 4 15 11 All units of length must be in km and time must be in seconds 6 Auxiliary Functions 6 1 deriv mex deriv mex is a function that can be used to obtain the output from a ODE function directly The syntax of the function call is dy deriv_mex odefun t y0 extras The inputs odefun y0 and extras are the same as would be used when calling one of the integrators The difference is that only a single time is sent t for which the derivative of the state at y 0 is computed dy is a column vector containing that derivative In MATLAB type help deriv_mex for information 6 2 UTC2JED m The function UTC2JED m can be used to convert from a UTC calendar date to a Julian Ephemeris Date JED JED is similar to Julian Date except it specifically is computed for JPL Coordinate Time which is similar to Dynamic Barycentric Time TDB This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help UTC2JED for information 6 3 JD2JED m The function JD2JED m can be used to convert from a Julian Date computed from UTC to a Julian Ephemeris Date JED This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help JD2JED for information 6 4 JED2UTC m The function JED2UTC m can be used to convert from a Julian E
14. Functions The mathematical models and requirements used in the ODE functions along with the abbreviations used in the function names are described in detail in the sections that follow 5 1 4 Two body Problem TwoBody The two body problem where a spacecraft orbits a single massive central body The acceleration due to the central body is di pr T2body pe 10 where r is the position vector of the spacecraft with respect to the planet center of mass r is the magnitude of that vector 5 1 2 Circular Restricted Three body Problem CRTBP In the circular restricted three body model there are two massive bodies in orbit about their mutual barycenter as depicted in Figure 1 To simplify the model each body orbits the barycenter in the same plane in perfectly circular orbits A spacecraft with infinitesimal mass experiences forces due to the gravitational influence of both bodies simultaneously and the two bodies are approximated as point masses The more massive body is labeled P and the other is P The coordinate frame has its origin at the barycenter and rotates with the two bodies so that P and P are always on the x axis with the positive x direction going from P to P2 The positive y axis is parallel to the velocity vector of P The three body gravitational parameter is called u ma k m ma SN where m is the mass of P and m is the mass of P5 The mass of the primaries is nondimensionalized so that the
15. When Solar Radiation Pressure is modeled TurboProp treats it as simply an adjustment to the Sun s gravitational parameter so the partials of the acceleration due to SRP with respect to the position of the spacecraft are lumped into the partials due to the Sun s third body gravity The partials with respect to the reflectance cr are explicity computed O fsnp A lsat PSR OcR m rosas 108 See Section 5 1 7 for a description of the notation 44 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 1 10 Two dimensional states _2D The model has only two dimensions of position and velocity 5 1 11 Unscented Kalman Filter Integrators uk In the Unscented Kalman Filter Julier and Uhlmann 1997 the unscented transformation generates 2L additional sigma points that must be propagated forward L is the number of filter estimated states To facilitate this process a single derivative function can be used to propagate all of the points in a single call to the integrator as opposed to 2L 1 calls These derivative functions are the same as their counterparts without the uk tag but they loop over multiple states in the yO array and compute derivatives for each Thus one call to the integration function will integrate all 2L 1 states 5 2 ODE Function Descriptions 5 2 1 TwoBody This ODE function is a model of the two body problem where a spacecraft orbits a single massive central body The state y0 is a 6 x
16. are changes made since v 3 3 TurboProp is now compatible with the Python scripting language given dependencies described later The spherical harmonics formulation has been changed from the classical form in Vallado and McClain 2001 to the cartesian representation described in Gottlieb 1993 Although slightly not as computationally efficient it removes the singularity at the poles See the grav section for more details Introduced a derivative function and the associated code to model the gravity field based on a collection of mass concentrations mascon Removed dynamic allocation from all spherical harmonic code to increase code reliability and reduce computation time While changing the spherical harmonics model and adding the mascon model the interfaces to the Gravity directory functions were streamlined Added two new integrators one 4th and one 6th order symplectic Runge Kutta Added derivative functions that allow compatibility with the unscented Kalman filter These derivative functions are added on an as needed basis by the developers A graphical user interface GUI was added to TurboProp to assist in learning the software Improved execution speed and added vector input output capabilities to e1orb and randv Removed calls to getpos and fsetpos in the JPL ephemeris software This improves portability and reliability of the software Added the FillExtras c and FillExtras h files to centralize pop
17. be useful for any derivative functions requiring knowledge of the initial state or the possible integration times After instantiation the integration is performed by T Y Cl tspan y0 This command may be repeatedly called for different t span and y0 values 4 11 odefun Input odefun is a string with the name of the ODE function or the derivative function As an example assignment in MATLAB you could type odefun CRTBP Section 5 explains the ODE func tions included in TurboProp Unfortunately users can not easily create their own ODE functions to use with the TurboProp integrators Advanced programmers would have to code up their own derivative functions in C and compile them with their integrators For more information please contact Brandon Jones 4 12 tspan Input tspan is a vector either a row or a column with values of time In Python it is an array of floats The first value in t span has to be the time associated with the initial state given in y0 States will be output in Y for all of the values in t span Derivative functions may have required units for t span so verify that you are using the correct time scale The DE type derivative functions now require that tspan be in seconds where it used to be days The fixed step integrators also require you to create tspan so that the integration time step will evenly divide each interval in t span except the last 14 4 13 y0 Input 4 INTEGRATORS 4 13 y0
18. central body The central body also has a spherical harmonic gravity field and an exponential atmosphere model for drag The state y 0 includes a state transition matrix and is a 56 x 1 vector The first 7 values are yO 1 m position y0 2 y position y0 3 z position y0 4 gx velocity y0 5 y velocity y0 6 z velocity y0 7 Cp coefficient of drag y0 8 56 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields 47 5 2 ODE Function Descriptions 5 ODE FUNCTIONS extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras degordstm Section 4 15 7 extras gravfile Section 4 15 8 extras Am Section 4 15 9 extras EOP 2 vector Section 4 15 11 extras atmos Section 4 15 12 extras mu should contain ju the gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 5 2 7 TwoBody grav_ukf TwoBody grav ukf is a model of the two body problem where a spacecraft orbits a single massive central body However it is designed for use in an Unscented Kalman Filter UKF Julier and Uhlmann 1997 In the UKF the unscented transformation generates 2L additional
19. extras degord 2 will be extras degord 1 if it is larger than extras degord 1 The values in ext ras degord cannot be larger than the gravity field you supply in the gravity file 4 15 7 extras degordstm extras degordstm is a vector with two doubles converted to integers in C that is used by the grav stm type ODE functions extras degordstm 1 sets the maximum degree of the spher ical harmonic gravity model used in the variational equations ext ras degord 2 sets the maxi mum order of the spherical harmonic gravity model used in the variational equations For example if you wish to use a 20 x 20 gravity field in the variational equations used to compute the state transition matrix type the following code in MATLAB extras degordstm 20 20 If you wish to neglect all the aspherical gravity effects in the state transition matrix and just include the point mass contributions use code like the following in MATLAB extras degordstm 0 0 17 4 15 extras Input 4 INTEGRATORS If the values in extras degordstm are greater than those in ext ras degord than the value of extras degord will be used instead of using extras degordstm Likewise the effective value of extras degordstm 2 will be extras degordstm 1 if it is larger than extras degordstm 1 4 15 8 extras gravfile For the grav type ODE functions the name of the gravity field spherical harmonic coefficients or point masses file has to be sent to the Turb
20. for 1 Jan 2007 16 38 35 UTC would be extras EOP 1 2454110 extras EOP 2 0 194215092592593 EOPfind m will automatically compute these times when it creates the extras EOP vector Oth erwise Terrestrial Time conversions can be performed using UTC2TT m or JED2TT m Since inte gration in the DE type derivative functions is performed using JED as the time scale the Orient function will automatically convert from JED to JDrr 19 4 15 extras Input 4 INTEGRATORS The polar motion parameters x and y are computed for time t in JDrr using a linear interpolation Ep t Epltres Eplt tres 3 Yp t Ypltres p t tres 4 Earth rotation parameters are computed for time t in JD using a linear interpolation wMi UTG e UFI 0TO wr DTOU eds 5 LOD t LOD trep ropi d 6 IERS precession and nutation correction parameters are computed for time t in JDrr using a linear interpolation d OAWisso t 0AViogoltres d OA 9go t tres 7 d OAeigso t OAeigso treg qe Aerosol tres 8 Using this linear interpolation the resulting EOP are only accurate for about one day so new EOP will need to be computed and the integration restarted about every 24 hours Section 5 1 6 describes these parameters and how they are used in detail 4 15 12 extras atmos extras atmos contains three parameters used to compute the atmospheric density for drag compu tations extr
21. frame is at the three body barycenter extras mu must be included and should contain u the three body system gravitational parameter The units of length time and mass must be the nondimensional units described in Section 5 1 2 5 2 10 CRTBP stm CRTBP stm is a model of the circular restricted three body problem in a rotating frame described in Section 5 1 2 The state y0 is a 42 x 1 vector The first 6 values are 49 5 2 ODE Function Descriptions 5 ODE FUNCTIONS yO 1 am position LU y0 2 y position LU y0 3 z position LU y0 4 x velocity LU TU y0 5 y velocity LU TU y0 6 z velocity LU TU y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is at the three body barycenter extras mu must be included and should contain y the three body system gravitational parameter The units of length time and mass must be the nondimensional units described in Section 5 1 2 5 2 11 CRTBP_2D CRTBP_2D is a two dimensional model of the circular restricted three body problem or the planar CRTBP in a rotating frame described in Section 5 1 2 The state y0 is a4 x 1 vector yO 1 am position LU y0 2 y position LU y0 3 z velocity LU TU y0 4 y velocity LU TU The origin of the coordinate frame is at the three body barycenter extras mu must be included and should contain y the three body system gravitational parameter The units of le
22. it is not extras planets 4 1 Marsis used 0 Mars is not extras planets 5 1 Jupiter is used 0 Jupiter is not extras planets 6 Saturn is used 0 Saturn is not extras planets 7 1 Uranus is used 0 Uranus is not extras planets 8 1 Neptune is used 0 Neptune is not extras planets 9 1 Pluto is used 0 Pluto is not extras planets 10 Any value will be automatically set to 0 extras planets 11 1 Sunis used O Sun is not extras planets 12 1 Earth is used 0 Earth is not extras planets 13 1 Moon is used 0 Moon is not The Earth Moon Barycenter models the combined mass of the Earth and Moon at a single point the Barycenter and should not be used when the Earth or the Moon are used separately The following example shows what users would type in MATLAB if they wished to include the Sun Venus Earth the Moon and Jupiter in the propagation extras planets 01 00100000 11 1 4 15 3 extras radius extras radius is a scalar double used in the CRTBP_grav type ODE functions It contains the radius of the secondary body in length units LU extras radius is optional for the DE_grav type ODE functions and can be used to assign a customized value for the central body s reference radius in km This is the radius used in the spherical harmonic equations but it is not used for eclipse calculations 4 15 4 extras center extras center is a scalar double converted to an integer in C that specifies
23. odefun is a string with the following ephemeris types Earth Earth stm Moon Moon stm The terms Earth or Moon specify the planet upon which the station or vehicle is located The term _stm specifies that a state transition matrix STM is desired in the output matrix Y tspan is a vector of times in seconds for which the inertial state of the surface object is de sired yO is the inertial state vector of the ground object at the time indicated in tspan 86400 extras reftime The units must be km and km s The MATLAB function SurferIC m can be used to convert latitude longitude height heading and speed information into an inertial vector to input as the y0 vector options is a scalar that is O when the surface object is a ground station and 1 when it is a vehicle extras is a structure used to pass extra parameters to the Surfer function For Surfer m the following fields are required extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras reftime JED Section 4 15 10 extras EOP 15 vector Optional when extras center 11 Section 4 15 11 As an example on how to used extras center if a station were on the surface of the Moon but it s inertial coordinates and STM were desired in an Earth centered frame GCRF set extras center 11 and odefun Moon_stm 63 6 18 TimeFrame m 6 AUXILIARY FUNCTIONS Surfer m works in the following way For the first value in tspan to the Earth orientation or lunar ori
24. org e The compiler used to compile you native version of Python Note SWIG is commonly distributed on Linux systems and the compiler is required to install Python If either dependency is not met the installation script will return an error During installation you may get warnings associated with the SWIG generated files _wrap c These are typical with some installations and have not been known to cause any problems with the generated Python libraries Additionally some warnings may be returned for files in the Gravity directory These can mostly be ignored If you have any questions e mail Brandon Jones After the software is installed the location of the TurboProp software must be added to the PYTHONPATH environment variable for future use in Python software Additionally the TurboProp directory itself must be specified in the user s environment in the TURBOPROP_DIR variable The exact values of these variables will be specified at the end of the installation script 8 3 PYTHON INTERFACE DEFINITION Unlike MATLAB Python does not provide a convenient which command for loading the path to the TurboProp data files However the TurboProp Python libraries provide the absolute path to the Data directory For example import os import TurboProp extras ephemfile os path join TurboProp DataDir DE 405 will insert the absolute location of the DE 405 file into the ext ras dictionary under th
25. plad ete 224244 X he a Be a Sea Sea Sen dead 16 LAS ezeras Ja US 2 445 464 664i 8 ESO RHEE ADEE e 16 amp l34 extras DODLEE Su e A Oe OR em ry vues 16 415 5 extras ephemfile 240 2 44929994 9 4 9 E VOU E Vox 17 amp l3h extracdegord AA axe E wR ek P eos 17 4 15 eutras dedgordeb t 6 6 2 bk RR Rx daa Sha EUR Rv 17 115 estraio Orari sce eh eae ee aa eed 18 BIS ETA E Lu RS A EA AA Po Ras XU uud he e ded 18 B islestrag teftti e cr basados ad 3 deg RARER dee deg 18 JIXDIbDeStradS BOE oscars aa aa a as oe ES OBS 18 a a es ux eR d v X V AOS BS Roe NUES 20 4 l3 l exztras attitude gt rrera Vps wo On em Re EO lee e da 21 a IB TETAS d Reve y 21 ZI GNI qoare qiix BOK Bo Gig Nod edic edic dor GE A RA 21 CONTENTS CONTENTS 5 ODE Functions Mathematical Theory for ODE Functions 22x yox RR Rs AE 5 2 Sd 5 1 2 53 13 5 1 4 24 3 2 8 9 3 2 1 8 dd 5 1 10 5 1 11 Two body Problem TwoBodw 222i 446 eh RR Romo Circular Restricted Three body Problem CRTBP JPL Ephemeris System DES escusa os dr AAA Mass Concentration Gravity Field PointMasses Spherical Harmonic Gravity Field grav o o Inertial Fixed Conversions 2 ee Solar Radiation Pressure SRP Atmospheric Drag drag caos ARA State Transition Matrix _stm ee ee ee Two dimensional states _2D eA Unscented Kalman Filter Integrators
26. point masses The state y0 is a 56 x 1 vector The first seven parameters are yO 1 x position km y position km z position km x velocity km s y velocity km s z velocity km s Cr reflectance y0 8 56 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras extras extras extras extras extras extras extras extras extras EOP 15 vector extras planets Section 4 15 2 Center Section 4 15 4 ephemfile Section 4 15 5 degord Section 4 15 6 degordstm Section 4 15 7 gravfile Section 4 15 8 A m m kg Section 4 15 9 reftime JED Section 4 15 10 mu km s optional Section 4 15 1 radius km optional Section 4 15 3 Optional when extras center 11 Section 4 15 11 56 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 2 23 DE drag grav DE drag grav is a model of the solar system using the JPL DE ephemerides atmospheric drag and a spherical harmonic gravity model for one of the planets The state y0 is a 7 x 1 vector yO 1 y0 7 x position km y position km
27. sigma points that must be propagated forward L is the number of filter estimated states To facilitate this process a single derivative function is used to propagate all of the points in a single call to the integrator as opposed to 2L 1 calls The central body also has a spherical harmonic gravity field The state y0 is a 78 x 1 vector yO 1 i 1 6 x position of the 7 th state y0 2 i 1 6 y position of the 7 th state y0 3 i 1 6 z position of the 7 th state y0 4 1 1 x6 x velocity of the th state y0 5 i 1 6 y velocity of the 7 th state y0 6 i 1 6 2 velocity of the i th state In this case 2 1 13 The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras EOP 2 vector Section 4 15 11 extras mu should contain ju the gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 5 2 8 TwoBody drag grav ukf TwoBody drag grav ukf is a model of the two body problem where a spacecraft orbits a single massive central body The central body also has a spherical harmonic gravity field and an exponential 48 5 2 ODE Function Desc
28. the ODE function the ODE function will ignore that field In Python extras is a dictionary with keys set to the fields of the extras structure For example extras gravfile os path join TurboProp DataDir GGMO2C sha Otherwise the definitions are the same For the following subsections this difference will not be repeated Any reference to extras x is equivalent to extras 4 15 1 extras mu extras mu is a scalar double containing a gravitational parameter For the CRTBP type ODE func tions this is the three body gravitational parameter or mass ratio For the TwoBody type ODE func tions this is the gravitational parameter of the central body extras mu is optional for the DE type ODE functions and can be used to assign a customized value for the central body s gravitational parameter in km s 15 4 15 extras Input 4 INTEGRATORS 4 15 2 extras planets extras planets isa vector of 13 values converted to integers used in the DE type ODE functions This vector indicates which planets will be included in the model A one indicates that the gravitational influence of a planet is used in the propagation and a zero or any other integer indicates that it is not Each parameter in ext ras planets looks like this extras planets 1 1 Mercury is used 0 Mercury is not extras planets 2 1 Venus is used 0 Venus is not extras planets 3 1 Earth Moon Barycenter is used 0
29. the potential in cartesian coordinates As derived in Gottlieb 1993 the following equations summarize the algorithm ay gt 5 gt mAnm Cus Eg T aim 26a n 2 m 1 5 x 3 i x m m SnmEm 1 an Cn mEm 1 26b n 2 m 1 a3 Y 3 Anal CnmEm F Bua 26c n 2 m 0 a4 Y 5 gt m DAnm CnmEm s Sam 26d n 2 m 0 u Z r a r p Zas ai 7 i 26e 28 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS where Dos Unm i Unm ye E ine m 2 dom 27 since n m um ls m Y 2 dom 2n 1 28 In the software I m is precomputed and stored at initialization The equations of motion are integrated in the inertial frame so the spacecraft state vector must be converted from the inertial frame to the body fixed frame before computing the accelerations After computing the body fixed accelerations they must be converted back into the inertial frame The conversion from an inertial to a fixed frame is described in detail for both the Moon and the Earth in Section 5 1 6 5 1 6 Inertial Fixed Conversions Lunar Orientation The conversion between inertial and fixed position for the Moon is computed using the Lunar Librations from the DE403 or DE405 ephemeris The DE403 ephemeris is best for the Moon and using Julian Ephemeris Date JED is technically more correct than using Terrestrial Time TT when interpolating the Lunar Librations The Lunar Librations are three Eu
30. to convert from an inertial coordinate system to an Earth or Planet fixed coordinate system For the TwoBody grav type functions ext ras EOP is a vector containing two doubles extras EOP 1 o Sidereal Time at time zero to extras EOP 2 Rotation rate of the central body These parameters are used to find the angle between the x axis of the inertial frame and the x axis of the surface fixed frame This angle is measured as a counter clockwise rotation about the positive z axis right hand rule the same as Greenwich Mean Sidereal Time The equation used to compute 0 t the sidereal time at time t is shown below A t 00 t 2 18 4 15 extras Input 4 INTEGRATORS Both quantities in extras EOP use units of radians to describe angles while ext ras EOP 2 uses whatever time unit is used in the state vector gravitational parameter and time span vector It is up to the user to make sure all the units are consistent For the DE_grav type functions where the central body is Earth ext ras EOP is a vector used to convert from an Earth centered Inertial frame GCRF to the Earth fixed frame ITRF This is done in the Orient function within the C code in the MEX functions and also in the EarthFrame m function in MATLAB They both use the 1976 Precession the 1980 Nutation with corrections Earth rotation and Polar Motion For details see Section 5 1 6 The extras EOP array can be filled auto matically using th
31. z velocity The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras EOP 2 vector Section 4 15 11 extras mu should contain ju the gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 5 2 4 TwoBody grav stm TwoBody grav stmis a model of the two body problem where a spacecraft orbits a single massive central body The central body also has a spherical harmonic gravity field The state y0 includes the state transition matrix and is a 42 x 1 vector The first 6 values are yO 1 ac position y0 2 y position y0 3 z position y0 4 az velocity y0 5 y velocity y0 6 z velocity y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras degordstm Section 4 15 7 extras gravfile Section 4 15 8 extras EOP 2 vector Section 4 15 11 extras mu should contain ju the
32. 0 GL derine uud unum a bus ar ecd ADU n Roe CUR un og ea ode Ee 60 C2 MUECA 6 keh ee ee oe a e eee 60 A ATP 60 hod JED UTICA AI 60 GS A II A 60 Bb UIC T a a a de e da eRe eS ere id ds 61 Ce Se E O D 61 5 9 JIEDETTAN occitano dd a i a de 61 CS GP ouvir A A a A eS we ud 61 BIO CRIBP DE atl lt i 0 d4 aa b medo de hee e o ERS GS 61 DII roti herte 2o AY OR A ORE OAR ORE SRR a REE Be 61 Gl TheseBo0 SUELE lt lt an Soe Aaa A Pe Ra ee OR S 61 D I3 BOPIVIQU B 2 heh he Ae See ee ee Peek eo Sk ede Rez ewe ede 62 6 14 WAGER rael tbh heh ek ROR aaa PHRASE WEES EE SES eS 62 6 15 Moon rame Me 242 243 SRE EAE RR ESRD ew eRe SOR Ewe me 62 OIO SUES Oa unde Bh oh Bd ke Sah Hoh ta aoe baad 62 MEN i er hf E Se ok ee oe ee ee GB ee BS Gn Ge ee eek 63 e suc dcr ee AGS SoG Oe Gok Gad Si bse Be old 64 A ee we ee eb Se a oa Bg ee w ex SS OS od owes 65 A 2 m DEC 65 7 MATLAB Graphical User Interface 65 1 TurboProp TurboProp is a package of functions that can be used to quickly propagate precise orbital or surface trajectories in the convenience of the MATLAB and Python programming environments To speed up the computations the orbit propagators are initial value ODE solvers coded in the C programming language The orbit propagators can be called as MATLAB functions through a MEX interface or Python via the SWIG module so the user gets both the superior execution speed of C code along with the ease of matrix computations and plotting that c
33. 0 ai ls Ox m Oys r 023 r3 r3 38 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS 23 Y3 and z3 are the inertial position vector from the center of mass of the third body to the space craft and r3 is the magnitude of that vector u3 is the gravitational parameter of the third body These computations are repeated for each third body included in the equations of motion Notice that this is the same as for the central body CRTBP STM For the Circular restricted Three body Problem the partials used to create matrix A are O CRTBP Ox OXCRTBP Oy OiCRTBP Oz OLCRTBP Oy OUCRTBP Ox OVCRTBP Oy O YorT BP Oz O VorTBP Ox OZcRTBP u li H 3 2n s OicRTBP Ox OZCRTBP Oz OUCRTBP Oy OZcRTBP Oz 3uz 31 2 H 1 p s Oz 3 r 5 2 5 3 T2 ri ri See Section 5 1 2 for details on the notation For the Circular restricted Three body Problem centered at the secondary with the 7 axis pointed toward the primary the partials used to create matrix A are 39 82 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS onrBP OL CRTBP Oy XonrBP OZ OTCRTBP Oy OYCRTBP OZ OYCRTBP Oy OYCRTBP OZ OVCRTBP m Od OZORTBP OZ OZCRTBP o OZCRTBP OZ 83 l p l p 1 p p i 3 3 5 5 mn ry ry T2 13 l p 1 pz 39 5 5 ri T3 3 1 4 1 u E 5 2
34. 1 Note the equation for S corrects a mistake in Gottlieb 1993 Finally Va Mn m nm E 4m l n 4 m4 n m n m 1 2 mo 100 Once the matrix is computed it must be rotated from the body fixed frame to the inertial frame r before being added to the rest of the A matrix This is done with the rotation y described earlier as 41 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS E J shown in Eqn 4 54 in Long et al 1989 The inertial version of is called ia the example Or OrGCRF below 3s m TGCRF ER y 101 OrGORF 3x3 Or 3 5 Ji CGCEF can then be added to the proper portion of the A matrix OVGCRF Drag STM For the TwoBody drag type derivative functions The variational equations used to create matrix A have contributions due to drag as shown below OF drug E Cp A Va x r Zm rH val O drag E Cp A mu ly Oy 2 m O drag u Cp A i Oz 2 mas rH OF aras ay Ve m ANA i has _ ET d PUT Oy 2 Vl OP drag Cp A pV Va z Oz 2 ara Obaras 1 A Cp 57 Va 2 Va O drag Cp A a ES Or m 2 5 2 V Val F Val OU drag e Cp A A ly Oy 2 m e Oljarag u Cp A e Oz 2 me w rH Olaas CoA Vas 1 8i 2 m 195 Oljarag e les ij 0j IVa Parag A 2 Vez V PE 2 m Va OU area u 1A C5 T 2m Vala 42 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS For the OT gu
35. 6 LAU Precession model given in Eqn 3 56 in Vallado and McClain 2001 C 2306 2181 Trr 0 301881 0 017998T 7 s 2004 3109 Trr 0 42665 T2 0 041833 T4 37 z 2306 2181 Tr7 1 09468 T7 0 018203T77 The units must be changed to radians Ty is the number of Julian centuries from J2000 in Terrestrial Time TT J Derr 2451545 0 i 36525 The equations for O and z are nominally written using 7 gp but Tr can be used with no discernible change in results The conversion for position and velocity from GCRF to MOD is Vallado and McClain 2001 Eqn 3 57 Trr 38 YmoD ROT3 z ROT 2 0 ROT3 G racnr VMOD ROT3 z ROT2 0 ROT3 G vaonr 39 Note that this MOD frame is different than the MOD frame obtained when converting from the Mean J2000 inertial frame The three rotation matrices shown above are combined and called PRE PRE ROT3 z ROT2 0 ROT3 C 40 1980 IAU Theory of Nutation corrections Next a conversion is performed from inertial MOD r yop to inertial True of Date TOD coordinates rrop This is done in TurboProp using the 1980 IAU Theory of Nutation First compute the mean obliquity Vallado and McClain 2001 Eqn 3 52 23 439291 0 0130042T pp 1 64 x 10 Tf 5 04 x 10 T3 41 The nutation in longitude AW and the nutation in obliquity Ae are normally computed from a se ries of 106 terms However TurboProp uses the interpol
36. 63 sin 2Q 49 AsT T OGMST F EQequinoz 50 Qq can be found using Eqn 3 54 in Vallado and McClain 2001 Q 125 04455501 1934 1361851T r 0020756174 4 2 139 x 107 Tp 1 650 x 10 Tf 51 The conversion from the TOD to a Pseudo Earth fixed PEF is Vallado and McClain 2001 Eqn 3 70 rper ROT3 9asr rrop 52 ROT3 04s7 will be called ST For the velocity the rotation rate of the Earth must be computed using the excess Length Of Day LOD Vallado and McClain 2001 Eqn 3 69 We 7 29211514670698 x 10 1 LOD 53 LOD is a parameter in the extras EOP vector see Section 4 15 11 and has units of days For units of seconds then the term in the equation above must be changed to LOD 86400 Using wa the conversion matrix ST is wesin Pasr wecos asr 0 ST we cos fasr wesin 0asr 0 54 0 0 0 Finally the velocity is transformed from TOD to PEF VpEr ST vrop ST rrop 55 Polar Motion To convert from PEF to the Earth fixed International Terrestrial Reference Frame ITRF a polar motion rotation must be applied x and y are the coordinates describing the motion of the Earth s rotation axis with respect to the Earth s crust These polar motion parameters are computed from the ext ras EOP vector see Section 4 15 11 Using x and yp the conversion from PEF to ITRF is Vallado and McClain 2001 Eqn 3 68 1 0 tp TITRE 0 1 w rrzr 56 XLp Yp 1 The matrix above wi
37. 64 6 19 randv m 7 MATLAB GRAPHICAL USER INTERFACE 6 19 randv m The MATLAB function randv m will convert Keplerian Orbit Elements to Cartesian coordinates For more information type help randv in MATLAB 6 20 elorb m The elorb m function will convert Cartesian coordinates to Keplerian Orbit elements For informa tion type help elorb in MATLAB 7 MATLAB Graphical User Interface In order to make learning TurboProp in MATLAB easier a GUI has been developed to automatically generate code This code is written to the MATLAB command window and can be copied into a separate file for later execution eo Figure 1 TurboProp v 3 3 SpoolUp Select your desired planetary model H Figure 3 Initial MATLAB GUI To execute the GUI run the command SpoolUp from the MATLAB command line The GUI should open and is depicted in Figure 3 The GUI is comprised of several elements that only activated when other elements allow their usage For example the Gravity Field pane will only be enabled when the primary body is the Earth or the Moon see Figure 4 To begin select a desired planetary model from the drop menu in the top left corner After the options have been selected for the propagated orbit click the Write Code button The generated code will be written to the MATLAB command window for future use 65 REFERENCES REFERENCES The MATLAB GUI does not support all of the capabilities of TurboProp For example the r
38. BP barycentric coordinate system to a CRTBP grav stm coordinate system centered on the secondary with the x axis positive in the direc tion of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius LU Section 4 15 3 extras degord Section 4 15 6 extras degordstm Section 4 15 7 extras gravfile Section 4 15 8 extras mu should contain u the three body system gravitational parameter see Section 4 15 1 extras radius should be included also and should contain the radius of the secondary body in length units LU see Section 4 15 3 All units of length time and mass must be the nondimensional units described in Section 5 1 2 5 2 15 DE DE is a model of the solar system using the JPL DE ephemerides The state y0 is a 6 x 1 vector yO 1 az position km y0 2 y position km y0 3 z position km y0 4 x velocity km s y0 5 y velocity km s y0 6 z velocity km s The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras reftime JED Section 4 15 10 ex
39. OTorTBpP Oy 1 7 1 hs De 5 Ti Ti T2 T5 l T s 22 m Ti T2 2 icnTBP _ Oz OUCRTBP Oz 1 1 Ez 1 5 ri ri T3 Pa See Section 5 1 2 for details on the notation Gravity Field STM To include the gravity field in the variational equations used to create the matrix A the partial of the gravity field portion of r with respect to r is computed These variational equations were derived in Gottlieb 1993 and are provided here as a reference The resulting matrix is Or pu 57s amp p e T pT PE bec Q kT E dl E v 40 N A Q Q N A R R A 84 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS where e 85 P k 0 0 1 E 86 d e Q R oj s T 0 87 A a3 a4 88 F L e Me 2 P a3 A 89 G Me P c a3 90 Additionally oo R n n L De 2 n m 1 n m 2 As Cnm Em Sum Fm 91 n 2 m 0 co R n n 7 n 2 m 0 oo R n TL 7 A 7 E E N gt 5 3 m m 1 Anm ComEm 2 Ss Fa 93 n 2 m 2 co R n n 7 7 7 E 7 n 2 m 2 RV P x 3n TU LIA m al at Om Em F Sod 95 n 2 m 0 oo R n n 7 7 7 E 7 Q zm y MA nm41ln m CnmEm 1 Boni 96 n 2 m 1 oo R n n 7 7 7 E E R 2 2 gt M Anm 1Tam SamEm 1 T Camina 97 n 2 m 1 oo R n n a 7 E E a es 3 n m 1 mol m Cn nEm 1 Duaci 98 n 2 m 1 co R n n 7 7 E 7 E Te 3 3n m 1 MAn m Sn mEm 1 BE Chumi m i 99 n 2 m
40. a da O drag Oy OL ares EN O drag r OU gra Oy O drag Oz OT asas OCp Dd Cp A Vale Vay a 3 me oe A Oird T Cp A Va ly _ Vas Oy 2 m eee Val Org B Cp A Valz Oz 2 mover rH O drag Cp A Vaz OQ 03 m Q9 OP dran z a il Vay Oy 2 Pei LM EY v 0 iv dog u d Ip Zm eVel A 500Cp pV m Voz W2 Vaz DE_drag type derivative functions The variational equations used to create matrix A have contributions due to drag as shown below rH Valy A 500C p p Van m rH Valz A 500 5 p View m A 500Cp p m A 500C p pV m A 500Cp pVaax m V2 ajo V va rH Va Va A 500C p pV V m 43 3 Va w Va od Va V w3 Vaz y Vaz a x UA Wy Vag ud Va Vos 105 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS thro 500C pp As mp cnl Vay EG ves Tes B 500C C pV EE eee hee 500C p p v ERE w4 Va Vag rte d ar 5000p Van 106 Dian 500C T Vay ai Tw 5000 p Va Va Mad 5000p Y A Vo eee E 50009 p A Vasa d Desa sucsA wv I Duces m E 500p PVs 107 E 800C pa Mires 50009 p Fes vall se 804 vil See Section 5 1 8 for a description of the notation SRP STM
41. alled DE 403 and the binary file for DE405 is called DE 405 The binary ephemeris files need to be installed on each machine as well since the binary data is stored differently on different computers The LP100K gravity field coefficients file is called jg1100k1 sha and the LP150Q coefficients file is called 391150q sha while GLGM2 sha contains the GLGM 2 gravity field The LP150Q file has the permanent tide correction already applied to Ja and C55 The Grace Gravity Model GGMO2C coefficients file also has the permanent tide corrections already applied and is called GGMO2C sha The JGM 3 gravity file is called JGM3 sha and the WGS 84 files is WGS84 sha IERS data from April 2007 is in the finals data file After the installation script completes you may run testephem m to make sure that the JPL ephemeris files were stored properly and can be read accurately test ephem m is in the Ephemeris directory and will make two plots of the errors If any of the blue error points are outside of the red dashed lines there may be a problem in your files The command window output will also tell you if there are any erroneous points The DE403 and DE405 errors may appear to be different magnitudes but that is just because the DE405 truth file testpo 405 does not have as many digits of precision as the DE403 truth file testpo 403 After verifying that the ephemeris files are working properly with testephem m you can erase all of the ASCII epheme
42. are the radii of the primary bodies Table 6 shows the values TurboProp uses for the planet radii Note that extras radius does not override these values for eclipse computations Most radius values were taken from Appendix D of Vallado and McClain 2001 but the radii for Mercury Venus and Mars were taken from the DE405 ephemeris Earth Moon represents the Earth Moon barycenter instead of an eclipsing body so it is given a radius of zero Table 6 PLANETARY RADII FOR ECLIPSE COMPUTATIONS Planet radius km Mercury 2439 76 Venus 6052 3 Earth Moon 0 Mars 3397 515 Jupiter 71 492 0 Saturn 60 268 0 Uranus 25 559 0 Neptune 24 764 0 Pluto 1151 0 Sun 696 000 0 Earth 6378 1363 Moon 1738 0 Units are km s 35 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS The relations below were used in determining 1 S2 g F 5 SS e S Rt O 6 O 68 26 0 S lt S SO This is based on the assumption that lt q gt s If that is false it only negatively affects the computation of when the Sun is partially eclipsed 5 1 8 Atmospheric Drag drag The acceleration due to drag for the TwoBody type derivative functions is computed using the equa tions Cp A Xdrag 75 p Val Val E Cp A Ydrag Pag Va 69 z Cp A Zdrag PVaslVal where C is the coefficient of drag which is part of the state vector is the area to mass ratio of the spacecraft specified by ext
43. as atmos 1 pp the atmospheric reference density or nominal density extras atmos 2 fro the atmospheric reference radius or base radius extras atmos 3 4H the atmospheric normalizing factor or scale height The density model used is an exponential atmosphere model and the density p at radius r from the center of the body is plr mer 9 For the TwoBody type derivative functions the units in po ro and H must agree with the units used in the integration and extras mu extras A_m For example if the state vector had units of km and km s ext ras mu had units of km s and extras A_m had units of km kg the units for po should be kg km and the units for ry and H should be km For the DE type derivative functions the units for py must be kg m the units for ry must be km and the units for H must be km Sample values for these parameters can be found in Vallado and McClain 2001 Table 8 4 p 534 2nd ed Ist Prnt 20 4 16 T Output 5 ODE FUNCTIONS 4 15 13 extras attitude extras attitude contains three or four parameters used to represent the vehicle attitude This field has only been used for select applications so their use has not been fully defined The initialization software allows the extras attitude vector to be of length 3 or 4 to accommodate both Euler angles or modified Rodriguez parameters or quaternions 4 16 T Output T is a vector of times an array in Python which for a
44. ated nutation parameters from the JPL DE ephemerides instead of computing all 106 terms in the trigonometric series Due to errors in 1976 IAU Precession and 1980 IAU Nutation the International Earth Rotation Service IERS publishes cor rections to the nutation parameters called OAW 19g and OAe19so These corrections also account for 31 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS differences between the GCRF frame used in TurboProp and the mean J2000 frame normally used with the 1976 IAU Precession and 1980 IAU Nutation models These nutation correction parameters are computed from the extras EOP vector see Section 4 15 11 AVi9g9 AV OAV 989 Aero Ae OAE i980 42 Ae10980 The conversion from MOD to TOD is Yrop ROT1 e ROT3 AV 980 ROT 1 r mop VTOD ROT 1 ROT3 AW 99 ROT1 amp v uop 43 The rotation ROT 1 e ROT 3 AW 1980 ROT 1 is called NUT Earth Rotation or Sidereal Time This conversion accounts for the instantaneous orientation of the Earth about its rotation axis First Greenwich Mean Sidereal Time 061157 must be computed Vallado and McClain 2001 Eqn 3 44 0ausr cusron weprecUT1 44 The UT 1 in the equation is the time past midnight in UT1 OGmsron is the Greenwich Mean Sidereal Time at midnight UT1 0h or 00 00 00 0 UT1 on the date of interest OGmsTon is computed using the equation Vallado and McClain 2001 Eqn 3 43 Oc
45. body problem in a rotating frame de scribed in Section 5 1 2 However CRTBP_grav_stm also includes a perturbation to the three body problem due to a spherical harmonic gravity model for the smaller primary or the secondary body Cur rently the software includes the assumption that the secondary body is in phase lock with the primary body i e the rotation period is the same as the orbit period Thus this function is only applicable to the Earth Moon system where the lunar gravity field is modeled The state y0 is a 42 x 1 vector The first 6 values are y0 1 g position LU y0 2 y position LU y0 3 position LU y0 4 amp velocity LU TU y0 5 y velocity LU TU y0 6 E velocity LU TU y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is not at the three body barycenter like the other CRTBP models but 51 5 2 ODE Function Descriptions 5 ODE FUNCTIONS the origin is at the center of the secondary with the x axis pointed toward the primary and the z axis perpendicular to the orbit plane of the primary and secondary Note that this reference frame is rotated 180 about the z axis from the other CRTBP models To convert from the standard CRTBP coordinates to the CRTBP_grav_stm coordinates use code like this y0 1 y0 1 1 mu y0 1 2 y0 1 2 5 y0 4 5 yO 4s5 7 This code converts a state vector yO from the normal CRT
46. e EOP find m function and a data file from IERS The file included with TurboProp finals data can be updated by downloading a new one from any of these locations ftp maia usno navy mil ser7 finals data ftp maia usno navy mil ser7 finals all http maia usno navy mil ser7 finals data http maia usno navy mil ser7 finals all For these Earth centered integrations in the JPL ephemeris system extras EOP is a vector con taining fifteen doubles extras HOP extras HOP 1 Integer portion of the EOP reference time trep in JDpr 2 Fractional portion of the EOP reference time tref in JD extras EOP 3 Nj Number of leap seconds at EOP reference time extras EOP 4 amp p tref bias term radians extras EOP 5 p rate term radians day extras EOP 6 Yp tref bias term radians extras EOP 7 Qp rate term radians day extras EOP 8 UT1 UTC tref bias term seconds extras EOP 9 UT1 UTC rate term seconds day extras EOP 10 LOD t ef bias term days extras EOP 11 LOD rate term days day extras EOP 12 0AViosoltre bias term radians extras EOP 13 LOAW 1980 rate term radians day extras EOP 14 OA ejog0 trer bias term radians extras EOP 15 ONe1980 rate term radians day To allow for sufficient numerical precision in some applications the reference time was split into two variables Thus the Julian Date in Terrestrial Time JD
47. e ephemfile key 3 Python Interface Definition Before moving on to the description of the integrators it is prudent to define the TurboProp interface when using the Python language If the user only plans on using the MATLAB capabilities skip to Section 4 This section assumes some familiarity with the Python language It is not intended to provide a tutorial or describe the language The user is directed to the many references on the web starting with the Python website http www python org TurboProp in Python was designed to exploit the object oriented capabilities and the modular defini tion of Python packages Once the software has been installed and the environment variables are set the basic TurboProp interface is provided by loading the TurboProp package i e import TurboProp This loads some basic definitions to make the interface with Python easier Specifically two variables are included in this top level package Dir andDataDir The Dir variable includes the absolute path to the TurboProp directory while the DataDir includes the absolute path to the Data directory Loading the main TurboProp module does not load the complete package which would require more time at initialization and loads potentially unnecessary software for a given application There are two packages within TurboProp that may prove useful to the user and they are accessed by import TurboProp Integrators import TurboPro
48. e vector in Y rk4fix dstm will reinitialize the state transition matrix STM portion of the state vector to an identity matrix Otherwise it is the same as rk4f ix 4 3 rk45 rk45 is a variable step Runge Kutta integrator that compares the results of 4th and 5th order integra tions to adjust the step size Like the ode45 function in MATLAB the integration scheme described in Dormand and Prince 1980 is utilized However the proprietary interpolation routine used by MAT LAB to determine the state at user specified nodes is not available so results may differ slightly when compared to ode45 rk45 is a variable step integrator meaning that the integration step size is ad justed to keep the difference between the fourth and fifth order Runge Kutta integrators less than a certain tolerance However the user still must input an initial guess or starting step size that will be used to begin the integration rk 45 requires that the opt ions vector have two values opt ions 1 contains the initial guess at the time step size or the step size used in the integration options 2 contains a tolerance used when comparing the seventh and eighth order Runge Kutta integrations The smaller this tolerance is the smaller the time steps must be so that the fourth and fifth order Runge Kutta integrations agree more closely to each other When TurboProp is compiled in debug mode a third option may be provided options 3 that determines if a file will be ge
49. emeris file for the positions and velocities of a planet Remember to input JED not the Julian Date based on UTC In MATLAB type help JPLDE for information 6 10 CRTBP DE m CRTBP DE m is used to convert an orbit generated in the CRTBP into the JPL ephemeris model centered on the Solar System barycenter or vice versa ThreeBodySystem m was changed to output LU and TU which can then be input directly into CRTBP DE m and rot inert m Now ThreeBodySystem m also offers automatic unit conversions In MATLAB type help CRTBP DE for information 6 11 rot inert m rot inert mis used to convert an orbit in the CRTBP into an inertial reference frame non rotating for the three body problem or vice versa The inputs to the function rot inert m were changed to be compatible with the new ThreeBodySystem m In MATLAB type help rot inert for information 6 12 ThreeBodySystem m ThreeBodySystem m is a function that can be used to obtain two and three body gravitational parameters the nondimensional length and time units and radii of the primaries for several three body systems ThreeBodySystem m was changed to output LU and TU which can then be input 61 6 13 EOPfind m 6 AUXILIARY FUNCTIONS directly into CRTBP DE m and rot_inert m In MATLAB type help ThreeBodySystem for information 6 13 EOPfind m The function EOP find m can be used to search an IERS data file to find the proper Earth Orientation Parameters to fil
50. entation conversions y to and 7 to are computed using the techniques shown in Section 5 1 6 The inertial state yO will be called Xjnertiai and it is converted to a planet fixed state x fixea for a moving ground vehicle like this vto O 3x3 o Alto to Xinertial to 1 1 3 0 3x3 is a 3 x 3 null matrix For each successive time t in t span the new fixed state is computed by assuming that the planet fixed velocity remains constant Z 3x3 t to Ulaxs 0 a RC X fizea to 114 X fixed t I 3 3 is a 3 x 3 identity matrix Then the inertial state at time t is Or Osx Xinerna lt Ey y t 7 Putting it all together gives WOT Olaxs Ere Ec ls Olives HO pOr Olas Hlsxs y to to which can be abbreviated by combining the three matrices y b Xinertial to 116 Xinertial t t to Xinertiat to 1 17 t to is the state transition matrix OT Oax Eos E al br M 118 HOP pO LOs Hlsxs V to to For a motionless station the equations are altered so that the planet fixed velocity is always zero The state transition matrix would be b t to 119 olt to KAT 0 3x3 to O 3x3 HO O sxs Ozxg O sxs and is used to compute the inertial state vector as shown above 6 18 TimeFrame m The TimeFrame m function is used to convert between different time scales such as UTC to Julian Date In MATLAB type help Timeframe for information
51. equires some extra computations but those computations are per formed only once per function call Note the normal integrators can be used with the proper derivative function to generate the STM for batch processors like the least squares filter While the ODE solvers are used to generate trajectories for orbiting spacecraft the MATLAB func tion Surfer mis used to propagate trajectories in inertial space for ground stations or vehicles SurferIC m is useful for generating initial conditions for Surfer m using inputs of latitude lon gitude height heading and speed For information on Surfer m see Section 6 This capability is not yet implemented in Python 10 41 rk4fix 4 INTEGRATORS 4 31 rk4fix rk4fix is a fixed step fourth order Runge Kutta integrator Since rk4f ix is a fixed step integrator the intervals in t span must be evenly divisible by the time step The only exception is the last interval which may be any length of time rk 4 ix mex requires that the options vector have two values options 1 contains the time step size or the step size used in the integration options 2 contains a tolerance used to verify that the time step evenly divides the time span intervals except the last If your time step evenly divides the time span but numerical issues make it so that it appears it does not use a larger value for options 2 A good value is normally 1e 6 4 0 rk4fix dstm After storing the output of the stat
52. ere r is the position vector of the spacecraft with respect to the planet s center of mass r is the magnitude of that vector u is the gravitational parameter of the body The acceleration due to other planetary bodies or third bodies Tapley et al 2004 is A Y3 sat To 3 P3body H3body Re F 3 J gt 19 3 sat TO where 13554 is the gravitational parameter of the third body r3 sat is the vector from the third body to the satellite rg 3 is the vector from the central body to the third body and 73 sat and re are the vector magnitudes By default the gravitational parameters are taken from the JPL ephemeris as shown in Table 1 The user can provide ext ras mu as an optional input to the DE type ODE functions and its value will be used as the gravitational parameter for the body that is the center of integration The gravita tional parameters for the other planets will remain at their default values The value in extras mu will also be used instead of the gravitational parameter in a gravity file if the ODE function requires it 24 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Table 1 DEFAULT GRAVITATIONAL PARAMETERS FROM JPL EPHEMERIDES Planet DE403 DE405 Mercury 2 20320804864179 x 104 2 20320804864179 x 104 Venus 3 2485859882646 x 10 3 2485859882646 x 10 Earth Moon 4 03503234715983x10 4 03503233479087 x 10 Mars 4 28283142580671x10 4 28283142580671 x10 Jupiter 1
53. eter Description Format Units Reference radius of the planet Scientific notation km Gravitational parameter Scientific notation km s Gravitational parameter uncertainty Degree of model field Order of model field Normalization Reference longitude normally 0 Reference latitude normally 0 Scientific notation Integer Integer Integer Scientific notation Scientific notation km s discarded N A discarded N A discarded 0 unnormalized 2 normalized degrees discarded degrees discarded Table 5 DATA LINE FORMAT FOR GRAVITY FIELD FILES C format string Sd Sd 1e 1e Sle 1e Parameter Description Format Degree index n Integer Order index m Integer Coefficient C Coefficient 5 Uncertainty in coefficient Chm Uncertainty in coefficient 5 Scientific Notation Scientific Notation Scientific Notation Scientific Notation The Earth and the Moon are the only planets supported in DE_grav type derivative functions so extras center should only be 11 or 12 for those functions The gravity field is modeled using U an aspherical Gottlieb 1993 potential function excluding the two body term xx RI 2 U xz y z 2395 5 An BH Cum Em Ss Fn 22 where r ai yb 23a Ey Ba Pu Epa Bre 23b Fn Fiimat BiFma Fo 0 A 23c 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Note this form is defined in terms of the Ca
54. ficient y0 8 56 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also 57 5 2 ODE Function Descriptions 5 ODE FUNCTIONS The ext ras structure should have all of the following fields extras extras extras extras extras extras extras extras 5 225 DE drag SRP grav planets extras extras extras extras center ephemfile degord degordstm gravfile Am reftime atmos mu radius Section 4 15 2 Section 4 15 4 Section 4 15 5 Section 4 15 6 Section 4 15 7 Section 4 15 8 m kg Section 4 15 9 JED Section 4 15 10 kg m km km Section 4 15 12 km s optional Section 4 15 1 km optional Section 4 15 3 EOP 15 vector Optional when extras center 11 Section 4 15 11 DE drag SRP grav is a model of the solar system using the JPL DE ephemerides atmospheric drag solar radiation pressure and a spherical harmonic gravity model for one of the planets The state yO is a 8 x 1 vector yO y0 8 x position km y position km z position km x velocity km s y velocity km s z velocity km s Cp drag coefficient Cr reflectance The o
55. gence tolerance is set in options 2 Unlike the explicit Runge Kutta integrators previously described symplectic integrators preserve the Hamiltonian of the system Thus a small time step is not required to preserve accuracy Of course these time steps vary with the system propagated More information on symplectic integrators can also be found in Hairer et al 2002 4 8 rk6sym The rk6sym is a like the rk4sym integrator except it uses the 3 stage Gauss Legendre scheme in Leimkuhler and Reich 2004 Given the higher order solution a larger time step than the rk4sym without a loss of accuracy 4 9 MATLAB Integrator Call Syntax The syntax for calling the TurboProp integrators was designed to be similar to the ODE solvers in MATLAB such as ode45 One of the formats for calling ode45 is 12 4 10 Python Integrator Call Syntax 4 INTEGRATORS T Y ode45 odefun tspan yO options extras The syntax for the TurboProp integrators is k4fix_mex odefun tspan yO options extras k4fix_dstm mex odefun tspan yO options extras HHHHHHHAH KKK KK KK i r 1 rk45_mex odefun tspan yO options extras rk45 dstm mex odefun tspan yO options extras rk78_mex odefun tspan yO options extras rk78_dstm_mex odefun tspan yO options extras rk4sym_mex odefun tspan yO options extras T Y rk 6sym mex odefun tspan yO0 options extras The inputs and outputs for these integrators a
56. gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 46 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 2 5 TwoBody drag grav TwoBody_drag_grav is a model of the two body problem where a spacecraft orbits a single mas sive central body The central body also has a spherical harmonic gravity field and an exponential atmosphere model for drag The state y0 is a 7 x 1 vector yO 1 m position y0 2 y position y0 3 z position y0 4 az velocity y0 5 y velocity y0 6 z velocity y0 7 Cp coefficient of drag The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras Am Section 4 15 9 extras EOP 2 vector Section 4 15 11 extras atmos Section 4 15 12 extras mu should contain ju the gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 5 2 6 TwoBody drag grav stm TwoBody drag grav stm is a model of the two body problem where a spacecraft orbits a single massive
57. he Moon finding the Moon fixed position is much easier since eq is zero and Cq Sq Re Re is set to 1737 4 km to be consistent with lunar topographical data although this is different than the radius used in gravity field computations 62 6 17 Surfer m 6 AUXILIARY FUNCTIONS For vehicles the heading and speed are converted to vsgz a velocity vector in SEZ coordinates South East Zenith cos a vsEz v sin a 111 0 v is the speed of the vehicle converted from km hr to km s a is the heading of the vehicle converted from degrees to radians and is zero for due north Vsgz is then converted to an planet fixed velocity V firea using Eqn 3 26 in Vallado and McClain 2001 sin ga cos A sin A cos ga cos A Vfized Sinl ya sin A cos A cos Qga sin A vsgz 112 cos dga 0 sin Hya T fixed AN v fizeq are then converted to the inertial frame using EarthFrame morMoonFrame m as appropriate 6 17 Surfer m Surfer mis a MATLAB function set up just like the ODE solvers except that it creates an ephemeris for a ground station or a ground vehicle instead of an orbiting satellite The ephemeris is created in the inertial frame just like satellite ephemerides so the ground stations can be used in orbit determination software without any special acommodation Here is its function call T Y Surfer odefun tspan yO options extras The outputs T and Y are just like the outputs of the ODE solvers
58. ices y for the Moon is y ROT3 v ROT1 0 ROT3 Q 33 Multiplying the three rotation matrices together gives CCH SCOS CHS SACOCHE SSO Y S v C 6 C v C 6 5 6 SS CW C AC e CSE 69 S 0 5 0 5 80 C 0 C 8 where C and S represent sine and cosine functions respectively The GEONS Mathematical Specification by Long and Lee 2006 gives the technique used to find y in Eqn 3 2 91a although different angles are used here The time derivative of y is taken assuming do _ do _ p a 0 and 7 0 Long and Lee 2006 give 1 a value of 0 2299708331 radians per day The conversion from Moon Centered Moon Fixed MCMF to Moon Centered inertial based on ICRF would be TICRF Di rucur ly vucur h rucur 36 VICRF 30 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Earth Orientation The conversion from an Earth centered inertial frame GCRF to an Earth cen tered Earth fixed frame ITRF is performed using the 1976 IAU Precession 1980 IAU Nutation with IERS corrections Earth rotation parameters and Polar Motion The conversion for velocity is included in the discussion even though it is not used to compute acceleration due to the gravity field 1976 IAU Precession To perform the conversion from inertial to fixed coordinates for the Earth a position vector in the inertial frame racrp must first be rotated into a Mean of Date MOD inertial frame This is done in TurboProp using the 197
59. inate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras Am m kg Section 4 15 9 extras reftime JED Section 4 15 10 extras mu km s optional Section 4 15 1 5 2 18 DE SRP stm DE SRP stmis a model of the solar system using the JPL DE ephemerides and solar radiation pressure The state transition matrix is integrated as well and includes contributions due to solar radiation pressure and the planetary point masses The state yO is a 49 x 1 vector The first 7 parameters are 53 5 2 ODE Function Descriptions 5 ODE FUNCTIONS yO 1 az position km y0 2 y position km y0 3 z position km y0 4 az velocity km s y0 5 y velocity km s y0 6 2 velocity km s y0 7 cpr reflectance y0 8 56 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also
60. ion calls The generation of the file is determined by the input options This capability allows the user to profile the performance of the integrator and assists with debugging However this creates some computation overhead that is not desired for most applications Thus this capability is only enabled when the software is compiled in debug mode In order to switch between debug mode and the non debug mode the installation script must be rerun with the appropriate option selected The installation scripts create the following files and directories in the TurboProp install directory lt InstallDir gt Data Ephemeris gravity and EOP data files Derivatives C derivative functions Ephemeris JPL ephemeris utilities Gravity Gravity field utilities Help MATLAB help files Integrators Integrators deriv_mex code MATLAB MATLAB functions PyUtils Python functions findptrmaker m Function to find ODE functions HelpFileMaker m Function to create help files InstallScript m Script to install TurboProp InstallScript py Script to install TurboProp TurboPath m Script to add TurboProp to the MATLAB path TurboTypes h C variable type definitions file init__ py Definition of the TurboProp package in Python The install scripts will also convert the ASCII ephemeris files into a more efficient set of binary 2 1 MATLAB Installation 2 INSTALLATION ephemeris files The binary ephemeris file for DE403 is c
61. k 4sym integrator is not listed as an option Code generated from SpoolUp can be adjusted accordingly to support all of the TurboProp features Future versions will support these additional features 80 Figure 1 TurboProp v 3 3 SpoolUp Select an htecrator htearator Options Single body Tnse Equator I Fixed step 4th order Runge Kutta Time Step 20 sec ASE O Variable 4 5th order Runge Kutta Time Step 106 STM Off m leranca z O Variable 7 8th order Runge Kutta _STMResot or C Derivative Only Initial State Primary Body s km Select the body forthe 6978 1363 d y 0 0 km Sun Y Planet Orientation 0 0 km Eann Theta rad rm 00 Time Span d 0 0 km sec Satum s Uranus terval 600 sec Neptune dTheta rad sec dy 6 64198160933358 km sec Pluto pe 7 29211585530 5 Time 3200 se az 3 60630177168135 km sec PE Gravity Field Spacecraft Properties CER ravity Field Off Drag Off nd SRP Off Figure 4 MATLAB GUI with the Gravity Field pane enabled Note Earth is selected as the primary body References Dormand J R and R J Prince A Family of Embedded Runge Kutta Formulae Journal of Compu tational and Applied Mathematics Volume 6 1980 pp 19 26 Gottlieb R G Fast Gravity Gravity Partials Normalized Gravity Gravity Gradient Torque and Magnetic Field Derivation Code and Data Technical Report NASA Contractor Report 188243 NASA Lyndon B Johnson Space Center Houst
62. k1 sha and the lunar LP150Q coefficients file is called j3g1150q sha jg1150q sha has the permanent tide correction already applied to J2 and C55 The GLGM 2 gravity field is from Lemoine et al 1997 and is called GLGM2 sha The gravity model for the Earth comes from the GRACE mission Tapley et al 2005 and is called GGMO2C sha The JGM 3 gravity field Tapley et al 1996 is called JGM3 sha while the WGS 84 gravity field NGA 2000 is named WGS84 sha The gravity field format was based after the file format for the lunar gravity field jg1100k1 sha found at http pds geosciences wustl edu missions lunarp shadr html The first line is a header and contains eight parameters three of which are used and the rest discarded The parameters are shown in Table 4 Note that these parameters are comma delimited If extras mu was not input by the user the gravitational parameter and reference radius from this header line will be used for all gravity calcula tions for the central body The user can provide extras mu and extras radius to override the parameters in the gravity header line After the header line data lines contain the gravity field coefficients These coefficients are comma delimited also and are described in Table 5 26 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Table 4 HEADER LINE FORMAT FOR GRAVITY FIELD FILES C format string Sle Sle Sle d d d Sle Sle Param
63. l harmonic gravity model for one of the planets The state yO is a 7 x 1 vector yO 1 az position km y0 2 y position km y0 3 z position km y0 4 x velocity km s y0 5 y velocity km s 6 z velocity km s y0 7 Cpr reflectance The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t sean must be in seconds also The ext ras structure should have all of the following fields 55 5 2 ODE Function Descriptions 5 ODE FUNCTIONS extras extras extras extras extras extras extras extras extras EOP 15 vector Optional when extras center 11 Section 4 15 11 extras planets Section 4 15 2 center Section 4 15 4 ephemfile Section 4 15 5 degord Section 4 15 6 gravfile Section 4 15 8 Am m kg Section 4 15 9 reftime JED Section 4 15 10 mu km s optional Section 4 15 1 radius km optional Section 4 15 3 5 2 22 DE SRP grav stm DE SRP grav stm is a model of the solar system using the JPL DE ephemerides solar radiation pressure and a spherical harmonic gravity model for one of the planets The state transition matrix is also integrated and the contributions of solar radiation pressure and a gravity model are included in the variational equations along with the planetary
64. l the extras EOP vector In MATLAB type help EOPfind for information 6 14 EarthFrame m The function EarthF rame m can be used with optional data from a ext ras EOP vector to convert an ephemeris between the coordinate systems GCRF J2000 MOD TOD and ITRF In MATLAB type help EarthFrame for information or see Section 5 1 6 6 15 MoonFrame m The function MoonF rame m can be used to convert an ephemeris between the Moon centered inertial ICRF and the Moon centered Moon fixed coordinate systems In MATLAB type help MoonFrame for information or see Section 5 1 6 6 16 SurferIC m The function SurferIC m can be used to generate initial conditions for Surfer m using inputs of latitude longitude height heading and speed In MATLAB type help SurferlIC for information For the Earth the Earth fixed position r7rrr is computed from geodetic latitude ga longitude A and ellipsoidal height herp using Eqn 3 14 in Vallado and McClain 2001 error in 2nd ed Ce hep co8 bga cos A TITRE Co hai cos gq sin A 109 Sa heup sin ga Cg and S are given from Eqn 3 7 in Vallado and McClain 2001 R Co E e sin ga 1 2 Se Re 1 e 110 y1 e2 sin boa Re is the equatorial radius of the Earth and is set to 6378 1363 km ea is the eccentricity of the Earth and Vallado and McClain 2001 gives it a value of 0 081819221456 on in section 3 2 of the 2nd ed For t
65. ler angles o 4 and 0 according to Newhall and Williams 1997 is the rotation angle along the Earth s mean equator of J2000 from the equinox to the ascending node of the lunar equator 0 is the inclination of the lunar equator with respect to the Earth s mean equator y is the angle along the lunar equator from the node to a longitude of zero Figure 2 from Newhall and Williams shows the angles graphically Make sure not to confuse the o mentioned here with the that represents latitude The conversion for velocity is included in the discussion even though it is not used to compute acceleration due to the gravity field Using these angles from the JPL ephemeris the conversion from the Moon Centered inertial position and velocity based on ICRF to Moon Centered Moon Fixed MCMF would be LCTMCMF YVICRF VMCMF YVICRF t YTICRF 29 The rotation matrices ROT 1 ROT 2 and ROT3 are 1 0 0 ROT1 0 0 cos 0 sin 0 30 0 sin 0 cos 0 cos 0 0 sin 0 ROT2 0 0 1 0 31 sin 9 0 cos 0 cos 0 sin 0 0 ROT3 0 sin 0 cos 0 O 32 0 0 1 29 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS nz 24 LI Lon J2000 Equinox Figure 2 The equatorial system in which the Libration Eular angles are defined The value of would be negative as shown credit Newhall and Williams 1997 Using these rotation matr
66. ll be called PM and the velocity conversion is virRF PM v per 57 33 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS To summarize the conversion from GCRF to ITRF ITRF PM ST NUT PREC rccorr Let y be PM ST NUT P REC and let y be PM ST NUT P REC and the notation is shortened and the conversion is YITRF recae Virrr lvracar y vacrr 59 Converting back to the inertial frame GCRF would be YGCRF Di rrrre vccnr WX rrrrr virre 60 5 1 7 Solar Radiation Pressure SRP A solar radiation pressure model is included in the force model for the satellite This is based on a constant area constant reflectance model The reflectance is also included in the spacecraft state vector The shadow model determines an approximate value of the percentage of the Sun s face visible at the spacecraft location The radii of the bodies used in the shadow model are taken from Vallado and McClain 2001 or the JPL ephemeris The value of solar radiation pressure is adjusted based on distance from the Sun The acceleration due to SRP is Tsnp ponen O Ont 61 m rosas Cr is the reflectivity of the spacecraft Ac is the cross sectional area of the spacecraft facing the Sun in square meters m is the mass of the spacecraft in kilograms r 5 is the vector from the center of the Sun to the spacecraft psg is the pressure of solar radiation in Pa This is about 4 53 x 10 Pa at one AU and is fo
67. ll the integrators in the current version of TurboProp should be identical to t span 4 17 Y Output Y is a matrix Each row in Y gives the value of the state vector in row form at the corresponding time in T For example Y 3 will have the state vector integrated to time T 3 which is the same as tspan 3 If the state transition matrix STM is part of the state vector it can be extracted in the following way when the state is a six vector STM reshape Y 1 6 1 6 36 6 6 In this example i is the row number The STM for time t span i would be stored in the variable STM as a6 x 6 matrix In Python Y is a list of lists where each list corresponds to a state at the corresponding time in the T array 5 ODE Functions Details on the state vectors and ext ras vectors for all the ODE functions are given later in this sec tion The current ODE functions in TurboProp are TwoBody CRTBP DE PointMasses TwoBody_stm CRTBP stm DE stm TwoBody_grav CRTBP_2D DE_SRP TwoBody grav stm CRTBP stm 2D DE SRP stm TwoBody_drag_grav CRTBP grav DE grav IwoBody drag grav stm CRTBP grav stm DE grav stm TwoBody grav ukf DE SRP grav TwoBody drag grav ukf DE SRP grav stm DE drag grav DE drag grav stm DE drag SRP grav DE drag SRP grav stm The list of accepted ODE functions can be obtained in MATLAB by typing help odefun 21 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS 5 1 Mathematical Theory for ODE
68. ly different For many applications only the function def inition is necessary However with each call to the TurboProp integration software some initialization is performed based on the odefcn used For example ephemeris files may be loaded into memory For repeated calls to the integrators this becomes inefficient The class interface is used to prevent this since these initialization tasks are only performed at the time of the class instantiation Assuming the same odefcn options and extras the integration function can be repeatedly called without 13 4 11 odefun Input 4 INTEGRATORS reinitializing the integrator Assuming the modules have been imported the calls to instantiate the classes are Cl rk4fix RungeKutta4Fix odefcn options extras tspan tspan yO y0 Cl rk4fix_dstm RungeKutta4Fix_dstm odefcn options extras tspan tspan yO y0 Cl rk45 RungeKutta45 odefcn options extras tspan tspan y0 y0 Cl rk45 dstm RungeKutta45 dstm odefcn options extras tspan tspan y0 y0 Cl rk78 RungeKutta78 odefcn options extras tspan tspan y0 y0 Cl rk78_dstm RungeKutta78_dstm odefcn options extras tspan tspan y0 y0 Cl rk4sym RungeKutta4Sym odefcn options extras tspan tspan yO y0 Cl rk6sym RungeKutta6Sym odefcn options extras tspan tspan yO y0 As indicated above congruent with all Python interface documentation the tspan and y0 inputs are optional at the time of class instantiation They may
69. mass of P is defined to be y and the mass of P is 1 u One nondimensional length unit LU is equal to the distance between the two primaries so the distance along the x axis from the origin to P is u LU and from the origin to P5 is 1 yy LU The time unit TU is defined such that P orbits around P in 27 TU The equations of motion for the CRTBP are Murray and Dermott 1999 z py gue 1 3 3 Ti T3 j 2 1 zh jv 12 ri T 1 zZ 5 ry T2 where ri zu y 2 and ro y z u 1 y 22 13 The functions CRTBP_DE m and rot_inert m can be used to convert from an inertial frame to the rotating non dimensional frame used in the CRTBP ODE functions See Section 6 for more details on those functions 2y 2 1 4 22 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Figure 1 Diagram of the Circular Restricted Three Body Problem with a rotating nondimensional coordinate frame The origin of the coordinate frame for the CRTBP_grav type functions is not at the three body barycenter like the other CRTBP models but the origin is at the center of the secondary with the x axis pointed toward the primary and the z axis perpendicular to the orbit plane of the primary and secondary Note that this reference frame is rotated 180 about the z axis from the other CRTBP models These new axes are called 7 y and 2 and the coordinates can be converted like this R
70. mston 100 4606184 36000 77005361TT1 0n 00038793177 on 2 6 x 10 Tri on 45 The units must be changed to radians Tiy71 0 is the number of Julian Centuries elapsed from J2000 to 00 00 00 0 UT1 on the date of interest J Duri on 2451545 0 Ti 46 UT1 0h 36595 46 J Dyin is the Julian Date at 00 00 00 0 UTI for the date of interest JDyr is computed from JD rr in the following way JDyri JDrr 32 184 UT1 UTC Mieap 86400 47 where 32 184 is the conversion from TT to TAI Nicap is the number of leap seconds needed to convert from TAI to UTC and UT1 UTC is the conversion from UTC to UTI Neap and UT1 UTC are computed from the EOP parameters in ext ras EOP see Section 4 15 11 Weprec 18 the Earth s mean angular rotation and can be computed using Eqn 3 38 in Vallado and McClain 2001 which gives revolutions per day Weprec 1 002737909350795 5 9006 x 107 Trr 5 9 x 107 PTE 48 32 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS The units must be converted to radians and 777 is the number of Julian centuries TT elapsed since J2000 The equation is nominally written using 77771 but Tr can be used with no discernible change in results After computing Oc sr it is necessary to compute the Greenwich Apparent Sidereal Time 0 Asr using the equation of the equinoxes Vallado and McClain 2001 Eqn 3 65 EQequinox AW cos 0 00264 sin Q 0 0000
71. nerated that includes each evaluation time of the derivative function and the total number of evaluations To enable this capability set options 3 equal to 1 1 0 in Python Otherwise a value of 0 or no value will pre vent the file from being generated The file will be created in the current working directory and named rk45out txt 4 4 rk45 dstm After storing the output of the state vector in Y rk 45 dstm will reinitialize the state transition matrix STM portion of the state vector to an identity matrix Otherwise it is the same as rk45 4 5 rk78 rk78 is a variable step Runge Kutta integrator that compares the results of 7th and 8th order inte grations to adjust the step size This integrator algorithm was modified by Jeff Parker using the code 11 46 rk78_dstm 4 INTEGRATORS ode78 v1 11 in MATLAB written by Marc Compere in 2001 rk78 is a variable step integra tor meaning that the integration step size is adjusted to keep the difference between the seventh and eighth order Runge Kutta integrators less than a certain tolerance However the user still must input an initial guess or starting step size that will be used to begin the integration rk78 requires that the options vector have two values options 1 contains the initial guess at the time step size or the step size used in the integration options 2 contains a tolerance used when comparing the seventh and eighth order Runge Kutta integrations The smaller this
72. ngth time and mass must be the nondimensional units described in Section 5 1 2 5 2 12 CRTBP stm 2D CRTBP stm 2D is a two dimensional model of the circular restricted three body problem or the pla nar CRTBP in a rotating frame described in Section 5 1 2 The state y0 is a 20 x 1 vector The first 4 values are yO 1 az position LU y0 2 y position LU y0 3 x velocity LU TU y0 4 y velocity LU TU y0 5 20 contain the state transition matrix in column order as described in Section 4 11 The origin of the coordinate frame is at the three body barycenter extras mu must be included and should contain y the three body system gravitational parameter The units of length time and mass must be the nondimensional units described in Section 5 1 2 5 2 13 CRTBP grav CRTBP grav is a model of the circular restricted three body problem in a rotating frame described in Section 5 1 2 However CRTBP grav also includes a perturbation to the three body problem due to a spherical harmonic gravity model for the smaller primary or the secondary body Currently the software includes the assumption that the secondary body is in phase lock with the primary body i e the rotation period is the same as the orbit period Thus this function is only applicable to the Earth Moon system where the lunar gravity field is modeled The state y0 is a 6 x 1 vector 50 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 1 g position LU 2 y po
73. oProp integrator This can be done in the following way extras gravfile which jgl150q sha 4 15 9 extras A m extras Amis a scalar double containing the area to mass ratio of the spacecraft The units depend on the ODE function used It is m kg for the DE type functions ext ras A_m is used in the solar radiation pressure SRP calculations The area is considered constant and is the cross sectional area of the spacecraft facing the Sun If ext ras A mis zero the SRP model is turned off 4 15 10 extras reftime For the DE type ODE functions the Julian Ephemeris Date JED is used to query the JPL DE eph emerides for planetary positions so the time scale of integration is in days To convert from a UTC calendar date to JED use the function UTC2 JED m described in Section 6 or JD2JED m Since JED numbers are so large some digits of precision are lost This can be remedied by sending a reference JED in the extras vector The time span in t span should be relative to that reference date In other words if the desired time span for the integration in JED was 2454003 5007 2454008 5007 The user would input a time span of tspan 86400x 0 5 and set the reference date extras reftime 2454003 5007 The scale factor of 86 400 is required because the time unit in DE type ODE functions is now seconds instead of days 4 15 11 extras EOP Earth Orientation Parameters EOP are contained in ext ras EOP This is used
74. of Geodesy Volume 79 8 2005 pp 467 478 Tapley B M Watkins J Ries G Davis R Eanes S Poole H Rim B Schutz C Shum R Nerem F Lerch J Marshall S Klosko N Pavlis and R Williamson The Joint Gravity Model 3 Journal of Geophysical Research Volume 101 B12 1996 pp 28 029 28 050 67 REFERENCES REFERENCES Tapley B D B E Schutz and G H Born Statistical Orbit Determination Elsevier Academic Press Burlington MA 2004 Vallado D and W McClain Fundamentals of Astrodynamics and Applications 2 Edition Micro cosm Press El Segundo CA 2001 68
75. ome with a high level interpreted language The orbit propagators or ODE solvers in TurboProp include a fixed step fourth order Runge Kutta solver a variable step fourth fifth order Runge Kutta solver a variable step seventh eighth order Runge Kutta solver and two symplectic Runge Kutta integrators of 4th and 6th order The surface sta tion vehicle propagator is a MATLAB function that generates the coordinates of the station vehicle in inertial space TurboProp was originally written to propagate orbits in the Circular Restricted Three body Problem and convert them to the JPL ephemeris Earth Orientation and a drag model have been added that will allow users to propagate orbits near Earth Included in this release are the JPL planetary ephemerides DE403 and DEA05 the lunar gravity models GLGM 2 LP100K and LP150Q and the Earth gravity models GGMO02C JGM 3 and WGS 84 1 TURBOPROP If you wish to keep the functions available for older versions of TurboProp please rename those functions so they won t conflict with the new ones in Version 4 0 Section 2 contains installation instructions for configuring TurboProp to work on your machine Sec tion 3 provides some details on the interface with Python Section 4 describes the integrators in detail Section 5 explains the ODE functions or derivative functions that are available with this version Section 6 explains how to use auxiliary MATLAB functions included with TurboProp The following
76. on TX February 1993 Hairer E C Lubich and G Wanner Geometric Numerical Integration Structure Preserving Algo rithms for Ordinary Differential Equations Number 31 in Springer Series in Computational Mathe matics Springer New York 2002 Hoffman D A A Set of C Utility Programs for Processing JPL Ephemeris Data NASA Johnson Space Center August 3 1998 Julier S J and J K Uhlmann A New Extension of the Kalman Filter to Nonlinear Systems Pro ceedings of SPIE Volume 3068 1997 pp 182 193 66 REFERENCES REFERENCES Konopliv A S W Asmar E Carranza W L Sjogren and D Yuan Recent Gravity Models as a Result of the Lunar Prospector Mission Icarus Volume 150 2001 pp 1 18 Leimkuhler B and S Reich Simulating Hamiltonian Dynamics 1 Edition Cambridge University Press New York 2004 Lemoine F G R D E Smith M T Zuber G A Neumann and D D Rowlands A 70th Degree Lunar Gravity Model GLGM 2 from Clementine and other tracking data Journal of Geophysical Research Volume 102 1997 pp 16 339 359 Long A C J O Cappellari C E Velez and A J Fuchs Goddard Trajectory Determination System GTDS Mathematical Theory Revision 1 July 1989 Long A C and T Lee GPS Enhanced Onboard Navigation System GEONS Mathematical Spec ifications Version 2 Release 2 6 Honeywell Technology Solutions Lanham MD February 23 2006 Lundberg
77. p PyUtils The Integrators package contains all of the integrators and the derivative function evaluator To access a specific integration module for example the rk4fix integrator use the import command import TurboProp Integrators rk4fix All of the rk4f ix integration capabilities will now be available in Python The same format is used for all integrators included in TurboProp For more information on the integrators and their interface 9 4 INTEGRATORS see Section 4 As dictated by the Python language the user can reassign the name of the imported integration module to another variable For example import TurboProp Integrators rk4fix as rk4fix Currently there are no modules of use for the user in the PyUt ils package Current modules are used by the software during installation or by the Integrators package The PyUtils package will include the Python equivalent to the MATLAB utilities described in Section 6 as they are developed and released Throughout this document Fortran indexing indices start with 1 is used to match the syntax of MATLAB However Python uses C indexing indices start with 0 To improve readability the MAT LAB notation is used throughout the document when talking about variables commonly used in both languages For Python users subtract the indices by 1 Additionally the user should make the appro priate changes to syntax y0 0 21 0 versus yO 1 1 0 4 Integrators
78. phemeris Date JED to a UTC calendar date This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help JED2UTC for information 6 5 JED2JD m The function JED2JD m can be used to convert from a Julian Ephemeris Date JED to a Julian Date based on UTC This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help JED2JD for information 60 6 6 UTC2TT m 6 AUXILIARY FUNCTIONS 6 6 UTC2TT m The function UTC2TT m can be used to convert from a UTC calendar date to a Julian Date in Terrestrial Time JDrr This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help UTC2TT for information 6 7 TT2JD m The function TT2JD m can be used to convert from a Julian Date based on Terrestrial Time JD to a Julian Date based on UTC This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help TT2JD for information 6 8 JED2TT m The function JED2TT m can be used to convert from a Julian Ephemeris Date JED to a Julian Date based on Terrestrial Time JD This function must be updated with the epochs of new leap seconds in Modified Julian Dates UTC In MATLAB type help JED2TT for information 6 9 JPLDE JPLDE is a function that can be used to query the JPL binary planetary eph
79. r Re Ma Nr Ho gi Re mne Il The converted equations of motion are 2 14 14 15 16 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS where ri vy 1 22 34 2 and ro i i 2 17 5 1 3 JPL Ephemeris System DE The JPL planetary ephemeris model DE ODE functions are in the ICRF inertial reference frame centered on a user specified planetary body or the Sun The reference plane of the ICRF is the mean equator of J2000 and the reference vector is the dynamic equinox of J2000 The positions and velocities of all the planets are obtained from JPL DE ephemerides DE403 Standish et al 1995 or DE405 Standish 1998 These ephemerides contain lunar orientation angles that can be used to rotate from the inertial Moon centered coordinates to Moon fixed coordinates The ephemerides also contain Earth Nutation angles for converting between inertial and fixed coordinates The binary ephemeris file for DE403 is called ephem 403 and the binary file for DE405 is called ephem 405 The code used to create the ephemeris files and to perform the interpolation was adapted from work by CCAR graduate student Greg Lehr and code by Mark Hoffman Hoffman 1998 at ftp ssd jpl nasa gov pub eph export C versions hoffman The ephemeris data files were obtained from ftp ssd jpl nasa gov pub eph export In the DE type functions the acceleration due to the central body is A r Y2body E 18 wh
80. ras A m V Vay and Vaz are the velocity components of the spacecraft with respect to the atmosphere assumed to rotate with the planet V is the magnitude of that vector The atmospheric density p at a certain radius from the center of the body r is computed using the parameters in extras atmos plr poe 70 See Section 4 15 12 for more information In TwoBody drag type derivative functions The atmospheric relative velocity vector is computed from the inertial velocities x y and 2 and positions x and y using the equations Vas iy Vy g x 71 Vaz is the rotation rate of the central body and is contained in ext ras EOP see Section 4 15 11 The acceleration due to drag for the DE type derivative functions is computed using the equations A Trad 500Cp mes Val A Udrag 500C p m P Ven Val 72 A Zdrag 7 500Cp mns Va g 36 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS The 500 coefficient is the conversion factor 1000 divided by two The conversion factor is needed because the units for 4 are m kg the units for p are kg m and the units for V are km s In DE_drag type derivative functions The atmospheric relative velocity vector is computed using the equation Vallado and McClain 2001 p 523 V T uxr 73 where r is the inertial velocity vector r is the inertial position vector and w is the Earth s rotation vector Wy 7 3 1 w
81. re described in detail in Sections 4 11 through 4 17 Note the calling function is the name of the integrator followed by the mex tag to prevent a conflict with any similarly named software on the MATLAB path 4 10 Python Integrator Call Syntax As previously mentioned the Python interface to TurboProp was designed to exploit the object oriented capabilities of Python To that end each TurboProp integrator has two access options a class interface and a function interface The function interface is intentionally similar to the MATLAB interface The integrator syntax using the shorthand defined in Section 3 is k4fix rk4fix odefun tspan y0 options extras k4fix_dstm rk4fix_dstm odefun tspan y0 options extras k45 rk45 odefun tspan yO options extras k45_dstm rk45_dstm odefun tspan yO options extras k78 rk78 odefun tspan yO options extras HAHAHAHAHA KKK KKK oH i k78_dstm rk78_dstm odefun tspan yO options extras k4sym rk4sym odefun tspan y0 options extras k6sym rk6sym odefun tspan y0 options extras T Y I I T I I E I E This syntax may seem redundant given the function name is the same as the module If the user wishes to prevent this redundancy the function can be assigned to a variable i e from TurboProp Integrators rk4fix import rk4fix as rk4fix T Y rk4fix odefun tspan yO options extras The class definition of the interface is slight
82. rigin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras extras extras extras extras extras extras planets center ephemfile degord extras extras extras extras mu gravfile Am reftime atmos radius EOP 15 vector Optional when extras center 11 Section 4 15 11 Section 4 15 2 Section 4 15 4 Section 4 15 5 Section 4 15 6 Section 4 15 8 m kg Section 4 15 9 JED Section 4 15 10 kg m km km Section 4 15 12 km s optional Section 4 15 1 km optional Section 4 15 3 58 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 5 2 26 DE drag SRP grav stm DE drag SRP grav stmis a model of the solar system using the JPL DE ephemerides atmospheric drag solar radiation pressure and a spherical harmonic gravity model for one of the planets The state transition matrix is also integrated and the contributions of a gravity model are included in the variational equations along with the planetary point masses solar radiation pressure and drag The state yO is a 72 x 1 vector The first eight parameters are i 1 c position km y position km 3 z position km 4 a velocity km s
83. riptions 5 ODE FUNCTIONS atmosphere model for drag This function is intended for use in a UKF The state y0 is a 105 x 1 vector yO 1 i 1 7 lt x position of the 7 th state y0 2 i 1 7 y position of the i th state y0 3 i 1 7 z position of the 7 th state yO 4 i 1 7 z velocity of the th state y0 5 1 1 x7 y velocity of the 7 th state y0 6 i 1 7 2 velocity of the 2 th state yO 7 i 1 7 Cp coefficient of drag of the i th state In this case 1 15 The origin of the coordinate frame is at the center of mass of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius Section 4 15 3 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras Am Section 4 15 9 extras EOP 2 vector Section 4 15 11 extras atmos Section 4 15 12 extras mu should contain ju the gravitational parameter of the central body see Section 4 15 1 extras radius should be included also and should contain the radius of the central body see Section 4 15 3 All units of length and time must be consistent 5 2 9 CRTBP CRTBP is a model of the circular restricted three body problem in a rotating frame described in Section 5 1 2 The state y0 is a 6 x 1 vector yO 1 am position LU y0 2 y position LU y0 3 z position LU y0 4 x velocity LU TU y0 5 y velocity LU TU y0 6 z velocity LU TU The origin of the coordinate
84. ris files to save hard drive space These files are named ascp1960 405 ascp1950 403 ascp1980 405 ascp1975 403 ascp2000 405 ascp2000 403 ascp2020 405 ascp2025 403 ascp2040 405 If you experience any difficulties or if you wish to have integration test scripts to verify that the inte grators are working properly contact Brandon Jones 2 1 MATLAB Installation In MATLAB run InstallScript m and it will compile the neccesary C code into MEX functions If this is your first time compiling a MEX function MATLAB may ask you which compiler you wish to use You can select the C compiler that comes with MATLAB which usually looks like 1 Lcc C version 2 4 Beyond the compiler used one can further customize the compilation of the C software by editing the mexopts sh file For example external libraries can be defined and included at compile time See the MATLAB documentation for more information When installing TurboProp on Linux boxes you may get a warning that the version of gcc on your machine is more recent than the versions used to test MEX compilation Usually this will not cause any problems but contact Brandon Jones if errors occur The new MEX functions will have a platform specific extension according to the table below 2 2 Python Installation 2 INSTALLATION System Type MEX File Extension Sun Solaris mexsol HP UX mexhpux Linux on PC mexglx Linux on AMD Opteron mexa64 Macintosh
85. rtesian not spherical coordinates Thus this formulation does not contain the singularity at the poles present in the classical formulation presented in Vallado and McClain 2001 and elsewhere R is the radius of the central body u is the gravitational parameter of the central body and n and m are the degree and order of the spherical harmonic model respectively Cash and Sra are normalized coefficients describing the magnitude of the spherical harmonics and are specific to the gravity field model being used A lo z is a normalized derived Legendre function computed using the following recursion Aoo Q 1 Aio Q av3 Ajo V3 1955227 Anna nm An 1n 1la n gt 1 2n Bon Anmlal BnmAn imlal An 2 m o m lt n n 1 m 24 where 2 1 2n 1 Ban c n Den 1 Um n m n m Since the B coefficients are not a function of a they are precomputed at software initialization and stored in memory for future use This formulation of the derived Legendre functions is adapted from personal correspondence with Dr Bob Schutz of the University of Texas Center for Space Research and from his book Satellite Dynamics Theory and Application current being written with G Giacaglia This adaptation to the formulation in Schutz s book was based on Lundberg and Schutz 1988 and is consistent with Gottlieb 1993 The acceleration in body fixed coordinates r is found by taking the gradient of
86. sition LU y0 3 Z position LU 4 amp velocity LU TU 5 y velocity LU TU y0 6 velocity LU TU The origin of the coordinate frame is not at the three body barycenter like the other CRTBP models but the origin is at the center of the secondary with the x axis pointed toward the primary and the z axis perpendicular to the orbit plane of the primary and secondary Note that this reference frame is rotated 180 about the z axis from the other CRTBP models To convert from the standard CRTBP coordinates to the CRTBP_grav coordinates use code like this y0 1 y0 1 1 mu y0 1 2 y0 1 2 y0 4 5 y0O 4 5 This code converts a state vector yO from the normal CRTBP barycentric coordinate system to a CRTBP grav coordinate system centered on the secondary with the x axis positive in the direction of the primary body The extras structure should have all of the following fields extras mu Section 4 15 1 extras radius LU Section 4 15 3 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras mu should contain u the three body system gravitational parameter see Section 4 15 1 extras radius should be included also and should contain the radius of the secondary body in length units LU see Section 4 15 3 All units of length time and mass must be the nondimensional units described in Section 5 1 2 5 20 44 CRTBP grav stm CRIBP_grav_stm is a model of the circular restricted three
87. system using the JPL DE ephemerides and a spherical harmonic gravity model for one of the planets The state transition matrix is also integrated and the contributions of a gravity model are included in the variational equations along with the planetary point masses The state y0 is a 42 x 1 vector The first six parameters are yO 1 az position km y0 2 y position km y0 3 z position km y0 4 az velocity km s y0 5 y velocity km s y0 6 2 velocity km s y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras degord Section 4 15 6 extras degordstm Section 4 15 7 extras gravfile Section 4 15 8 extras reftime JED Section 4 15 10 extras mu km s optional Section 4 15 1 extras radius km optional Section 4 15 3 extras EOP 15 vector Optional when extras center 11 Section 4 15 11 5 2 21 DE SRP grav DE SRP grav is a model of the solar system using the JPL DE ephemerides solar radiation pressure and a spherica
88. t is represented as a fraction of the total mass of the system The total mass is derived from the gravitation parameter given in the header line Locations of the point masses are defined in the planet fixed frame in spherical coordinates and rotated based on the ext ras EOP over time An example file example ptm is provided in the Data directory 25 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS Table 2 HEADER LINE FORMAT FOR POINT MASS FILES C format string Sle Sd Parameter Description Format Units Gravitation Parameter of the planet Scientific notation km s Number of point masses Integer N A Table 3 DATA LINE FORMAT FOR POINT MASS FILES C format string Sle 1e Sle 1e Parameter Description Format Units Fraction of total mass for given point Scientific Notation N A fraction of whole Point Radius Scientific Notation km Point Latitude Scientific Notation degrees Point Longitude Scientific Notation degrees 5 1 5 Spherical Harmonic Gravity Field grav A spherical harmonic gravity field model is used for one of the planetary bodies Current ODE func tions are only compatible with Earth and lunar gravity fields In other words Earth fixed and lunar fixed orientations are the only ones currently available Two lunar gravity fields are from Konopliv et al 2001 The lunar LP100K gravity field coefficients file is called g1100
89. the planet that is the center of the coordinate frame used in the integration according to the list below 16 4 15 extras Input 4 INTEGRATORS extras center 0 Mercury Centered Inertial extras center 1 Venus Centered Inertial extras center 2 Earth Moon Barycenter Centered Inertial extras center 3 Mars Centered Inertial extras center 4 Jupiter Centered Inertial extras center 5 Saturn Centered Inertial extras center 6 Uranus Centered Inertial extras center 7 Neptune Centered Inertial extras center 8 Pluto Centered Inertial extras center 10 Sun Centered Inertial extras center 11 Earth Centered Inertial extras center 12 Moon Centered Inertial extras center cannot be any other value 4 15 5 extras ephemfile extras ephemfile For the DE type ODE functions the name of the binary planetary ephemeris file has to be sent to the TurboProp integrator This can be done in the following way extras ephemfile which DE 405 4 15 6 extras degord extras degord is a vector with two doubles converted to integers in C that is used by the grav type ODE functions extras degord 1 is the maximum degree to be used in the spherical har monic gravity model extras degord 2 is the maximum order to be used in the spherical har monic gravity model For example if you wish to use a 50 x 50 gravity field type the following code in MATLAB extras degord 50 50 The effective value of
90. tolerance is the smaller the time steps must be so that the seventh and eighth order Runge Kutta integrations agree more closely to each other When TurboProp is compiled in debug mode a third option may be provided options 3 that determines if a file will be generated that includes each evaluation time of the derivative function and the total number of evaluations To enable this capability set options 3 equal to 1 1 0 in Python Otherwise a value of 0 or no value will prevent the file from being generated The file will be created in the current working directory and named rk780ut txt 4 6 rk78 dstm After storing the output of the state vector in Y xk 78 dstm will reinitialize the state transition matrix STM portion of the state vector to an identity matrix Otherwise it is the same as rk 78 4 7 rk4sym The rk4sym is a fixed step symplectic fourth order Runge Kutta integrator The integration scheme utilizes the 2 stage Gauss Legendre scheme describe in Leimkuhler and Reich 2004 Since rk4sym is a fixed step integrator the intervals in tspan must be evenly divisible by the time step The only exception is the last interval which may be any length of time rk4sym requires that the options vector have two values options 1 continas the time step size or the step size used in the inte gration Since the integration scheme is implicit a fixed point iteration is required to estimate the step using the time step The conver
91. tras mu km s optional Section 4 15 1 5 20 16 DE stm DE stm is a model of the solar system using the JPL DE ephemerides with a state transition matrix that includes the contributions of the planetary point masses The state y0 is a 42 x 1 vector The first 52 5 2 ODE Function Descriptions 5 ODE FUNCTIONS 6 values are yO 1 am position km y0 2 y position km y0 3 z position km y0 4 az velocity km s y0 5 y velocity km s y0 6 z velocity km s y0 7 42 contain the state transition matrix in column order as described in Section 4 11 The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t span must be in seconds also The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras reftime JED Section 4 15 10 extras mu km s optional Section 4 15 1 5 2 17 DE SRP DE SRP is a model of the solar system using the JPL DE ephemerides and solar radiation pressure The state y0 is a 7 x 1 vector yO 1 az position km y0 2 y position km y0 3 z position km y0 4 x velocity km s y0 5 y velocity km s y0 6 z velocity km s y0 7 cpr reflectance The origin of the inertial coord
92. ulation of the extras structure in MATLAB mex files Check for integrator failure after derivative function shutdown is complete Allow for options array of 3 values for debugging depends on compile mode and the integra tor used Add a check to see if the integration step size is a NaN or smaller than the machine precision Added an attitude field to the extras structure for some users 2 INSTALLATION e MATLAB time conversion software updated for the new leap second e Slight change to the eclipse fraction if test in SRP derivative functions This prevents the deriva tive function from generating a NaN when the satellite orbit radius is less than the radius of the primary body e Some general code clean up removal 2 Installation Before installing this new version of TurboProp you may wish to rename your older version TurboProp functions if you wish to retain the capability to use them You will need to download the zip file TurboProp4 0 zip containing TurboProp v 4 0 into the directory of your choosing Uncompress the zip file and run the follow the software specific installation instruction which follow in Sections 2 1 and 2 2 When running the installation script you will be prompted whether to compile in debug mode When the software is compiled in debug mode the rk45 and rk78 integrators can generate an output file that includes the time of each evaluation of the derivative function and the number of derivative funct
93. und by dividing the solar constant flux in W m by the speed of light p 544 Vallado and McClain 2001 1358 W m 62 PSRAU 299 792 458m s d For varying distances from the Sun psg would be 149 597 870 km DSR PSR AU l 63 rosal T Osa Should be in km A factor can be applied to represent the fraction of the Sun s face that is visible When km and km s are the units used for position and velocity respectively the acceleration due to SRP is 1 km TO sat 1000m rosal A orp cr 149 597 870km psg av 64 TL 34 5 1 Mathematical Theory for ODE Functions 5 ODE FUNCTIONS L 1 when there is a clear line of sight vector between the spacecraft and the Sun and 0 when the Sun s light is blocked by a body When determining the line of sight it was assumed that the center planet ext ras center is the only body that would eclipse the Sun The angle between the center of the planet and the center of the Sun is called and is found using the dot product m fg l sat 18 the vector from the center of the Sun to the spacecraft and rq sat is from the central body to the spacecraft sat and rq are the magnitudes of those vectors The angle between the center and limb for the two bodies are found using these equations assuming that the bodies are spherical CQOsat TC sat rosas Irq sat 65 R so sin O 66 nore R so sin K 67 Irq Re and Rg
94. z position km x velocity km s y velocity km s z velocity km s Cp drag coefficient The origin of the inertial coordinate frame is at the center of the specified celestial body using the ICRF coordinate system The unit of length must be kilometers and the unit of time is seconds so velocities are expressed in kilometers per second km s t sean must be in seconds also The ext ras structure should have all of the following fields extras planets Section 4 15 2 extras center Section 4 15 4 extras ephemfile Section 4 15 5 extras degord Section 4 15 6 extras gravfile Section 4 15 8 extras Am m kg Section 4 15 9 extras reftime JED Section 4 15 10 extras atmos kg m km km Section 4 15 12 extras mu km s optional Section 4 15 1 extras radius km optional Section 4 15 3 extras EOP 15 vector Optional when extras center 11 Section 4 15 11 5 224 DE drag grav stm DE drag grav stmis a model of the solar system using the JPL DE ephemerides atmospheric drag and a spherical harmonic gravity model for one of the planets The state transition matrix is also integrated and the contributions of a gravity model are included in the variational equations along with the planetary point masses and drag The state y0 is a 56 x 1 vector The first seven parameters are yO 1 x position km y position km z position km x velocity km s y velocity km s z velocity km s Cp drag coef

The User Manual can be found here. - CCAR

Contents

Download Pdf Manuals

Related Search

Related Contents