ABSTRACT Title of thesis: ORBIT DETERMINATION METHODS FOR DEEP SPACE DRAG-FREE CONTROLLED LASER INTERFEROMETRY MISSIONS Lauren Rose Chung, Master of Science, 2006 Thesis directed by: Assistant Professor Dr. Ella Atkins Department of Aerospace Engineering University of Maryland The objective of this thesis is to determine if different measurement combinations will improve absolute and relative orbit determination (OD) accuracy for the Laser Interferometer Space Antenna (LISA) mission. LISA, a triangular three-satellite heliocentric constellation, measures the relative length changes between spacecraft. The unique contribution of this work is the incorporation of inter-spacecraft range and range rate observations into the OD process, in addition to Deep Space Network (DSN) range and range rate, and Very Long Baseline Interferometry (VLBI) angle tracking data. MATLAB was used for orbit propagation, simulating measurements, and executing the weighted batch least squares algorithm to obtain the best state estimate. Covariance and Monte Carlo analyses were performed to compare the four test cases. The results concluded that cases with inter-spacecraft data had the least absolute and relative OD errors, VLBI eventually becomes ineffectual, and DSN needs at least 20 days of tracking to become observable. ORBIT DETERMINATION METHODS FOR DEEP SPACE DRAG-FREE CONTROLLED LASER INTERFEROMETRY MISSIONS By Lauren Rose Chung Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Masters of Science 2006 Advisory Committee: Assistant Professor Ella Atkins Dr. Tristram Tupper Hyde Associate Professor Robert Sanner ? Copyright by Lauren Rose Chung 2006 ii Preface This thesis topic began as a ten-week project during my Summer Institute in Engineering and Computer Applications (SIECA) internship at NASA-Goddard Space Flight Center during the summer of 2004. Under the mentorship of Dr. Tristram Tupper Hyde, the project involved taking a preliminary look at a two-dimensional model of a spacecraft around a flat Earth then simulating three spacecraft orbiting around the Sun in three- dimensions. The only data type considered was range. The batch least squares method was used to estimate the best state. The experience exposed me to brilliant people, got me involved in research, interested me in orbit determination, and made me want to keep on pursuing the complicated problem of tracking spacecraft. I returned to University of Maryland that fall to finish my Master?s coursework and in the spring of 2005, I fully embarked again on this journey of research to complete the final thesis requirement. The topic got deeper and many lessons were learned along the way. iii Dedication This thesis is dedicated to my Mother: Maxine "Joy" Chung who is simply mommy. This is a woman who is always so positive despite it all, who is so giving, who works so hard to provide the best for us, who has an endless supply of crazy expressions that taught me a lot about life, and most of all who always believed in me. iv Acknowledgements I first want to acknowledge all the members of my very large family, without the support of which, I could not have done this. The most special shout-out is to my sister Kristen who is my twin (okay, we are really one year and nine months apart) - all my quirkiness and humor is not lost on you. You get me and for that I am eternally grateful. Thank you Kristen for counseling me through all the stages of this journey: the pep talks and the straight-talks, and for keeping me hip and down-to-Earth when I thought I was losing my mind. Mom and Dad, I want to thank you for all the wake-up calls. Mom, thanks for always alerting me whenever the words space, shuttle, launch or NASA were mentioned on TV. You are that constant force behind me, for better or worse. Dad, thanks for loving science fiction even if I cannot get down with Stargate, Star Trek and Babylon 5. I know: "what kind of aerospace engineer am I?" However we do see eye to eye on our love of the Lord of the Rings and curiously Sex and the City, and thanks for helping me handle those other areas of life: finances and car issues. Thanks to my brother Hugh (an engineer in disguise who also loves Stargate and all those build-it shows) and his wife Shanna for being our get-away out in California. And to my bite-size inspirations, my niece and nephew: Mattea and Austin: the cutest kids in the world, who are so smart, funny and imaginative and who say the darnedest things. To my local family: my cousin Fran (and her dog Quincy) for helping me get acclimated to MD and Paul, Rona and Aunt Patsy (and Sam the cat) for being my escape in VA and v for listening while I went on and on about school, oh the drama. And thanks for feeding me Dave?s Barbeque. I want to acknowledge my friends for being there for me, supporting me, and somehow always knowing that I could do this. Special shout-outs to my most intelligent and beautiful friends: Keturah, Heneryatta and Angel who are doing their thing by getting their PhDs. You go girls!!! I cherish that we met at Penn State, formed a lifelong bond and had so many fun times in the process. You women continue to inspire me. Angel our friendship deepened here together at Maryland as we tried to navigate through this crazy stage called "Grad school life," and thanks for formally introducing me to "PhD Comics"- how relatable! More Penn Staters that I want to acknowledge are my dear friends Ebony and Brandon (and now baby Jada), Fernando- my ex-ballroom dance partner, Yanira, and Karla. Also, Amy Freeman, Sandra Johnson and Anita Persaud from the MEP, Diane and D?Andre, the II Hype Penn State NSBE Chapter, and my fellow Aerospace Engineers: Jonathan Coleman, Eric Parsons, Erika Mendoza and Dr. Melton for being such a fabulous teacher. Jen DeLecce, my partner in crime in Australia, who is only in VA, but might as well be Australia, since I barely make it to see you. Here is to more free time when this thesis is complete! And to my pal, Damon Bradley, our paths have coincided way too long; Penn State, NSBE ("II Hype!"), NASA internships, Grad school at UMD, and now living and working in Maryland. Damon thank you for always being so positive, smart, for converting me to love glow sticks and techno music, and for deriving crazy equations just for fun and making it okay to be a little (or a lot) geeky. vi I also want to acknowledge fellow UMD grad students, including Janisa Henry and my SSL lab mates, who are going through process, as well. And the members of Black Graduate Student Association (BGSA), especially for all the potlucks and deep conversations that went on for hours. Good luck to all. I also want to thank Cynthia Simmons for befriending me from the start in Astrodynamics and moving me to achieve. And thank you Dan Rolf for lending me (err? letting me kidnap) that Vallado textbook. I want to acknowledge the people in my hobbies/passions: my dance sisters (Malika, Sirena, Leyla, Zuhairah and Gia) and yoga buddies for adding a dose of drama to my life (the former) and for centering me (latter). I am truly thankful and grateful for the journey. I have learned so much more about myself and plan to dance and do yoga until I die. To all the mentors throughout my life who have inspired me every step of the way: Willard Manning who I credit with being my first true Mentor with a project called "The Environmental Benefits of Earth Observing Satellites," and who is now a life-long friend and unofficial family member. And to Sharae McKay who inspired me to go for a higher degree in Aerospace engineering. I would like to acknowledge these people from my NASA-GSFC 2004 internship: -Dr. Russell Carpenter: for answering all my questions, being so generous with your vast knowledge and being an inspiration. vii -Dr. Landis Markley: for letting me co-inhabit your office, being the human encyclopedia and pulling any journal article or resource I could possibly need from his files. From UMD: -Dr. Healy: though I met you late in the process, I thank you for your help this past summer, and your interest in Cathy Thornton' paper and class notes. -Dr. Bowden: for always being so positive every time I saw you and for your sense of amazement at everything your mere students were doing. -Dr. Tolson: for drilling into your student?s heads that a solar day is longer than a sidereal day, for assigning those group projects (my group: Dan, Matt and Adam) to make all those MATLAB Astronautics toolkits (though we students were not happy at the time) and for giving us that Cathy Thornton paper. -Dr. Atkins: for giving me a space in the Space Systems Lab, though I was the odd woman out in the robotics lab working on an orbit determination project! Thanks for your help and counseling on various stages of my thesis. And thank you immensely for your wonderful revisions of the thesis. -Dr. Sanner: for early conversations on orbit determination. And I would also like to thank the wonderful employees of a.i.-solutions, Inc. for hiring a girl, not quite done with her Master's thesis. Especially to David Hooshmand, Christy Jeffries, Conrad Schiff, Daryl Carrington who hired me. I won't disappoint you, I am just warming up. And my entertaining officemates: Allan, Dominic and Matt- we do have the hottest (temperature-wise) office at a.i. viii And once again, I want to sincerely thank my mentor Dr. Tristram "Tupper" Hyde who went above the call of duty and took on the task mentoring a student through her Master's thesis. You are the one who gave a hopeless grad student a Master's thesis topic in the summer of 2004 and stood by me for the last two years following that summer internship. Thank you for your genuine interest in this research topic. Thanks for the e-mails, the calls, the meetings, helping with the MATLAB codes in all its states of "bugginess," dealing with my endless reading, note taking and "research" that had to finally be channeled into concrete results. ix Table of Contents Preface........................................................................................................................... ii Dedication....................................................................................................................iii Acknowledgements...................................................................................................... iv Table of Contents......................................................................................................... ix List of Tables ............................................................................................................... xi List of Figures............................................................................................................. xii List of Abbreviations .................................................................................................. xv List of Symbols.......................................................................................................... xix Chapter 1: Introduction................................................................................................. 1 1.1. Thesis Objectives & Contributions.................................................................. 4 1.2. Overview of Thesis........................................................................................... 5 Chapter 2: Background ................................................................................................. 6 2.1. Mission & Science Objectives.......................................................................... 7 2.1.1. Gravitational Waves .................................................................................. 8 2.1.2. Laser Interferometry ............................................................................... 10 2.2. LISA Spacecraft Design ................................................................................. 11 2.3. Instrument Package ........................................................................................ 13 2.3.1. Drag-Free Control ................................................................................... 15 2.4. LISA Launch & Mission Evolution ............................................................... 16 2.5. Orbit Geometry & Configuration ................................................................... 17 2.6. Constraints & Technological Challenges........................................................ 20 2.7. Related Missions............................................................................................. 21 Chapter 3: Methods & Math Techniques.................................................................... 26 3.1. The Two-Body Problem ................................................................................. 26 3.2. Orbital Elements ............................................................................................. 29 3.3. Reference Frames............................................................................................ 34 3.3.1. Sun-Centered Inertial Coordinate System ............................................... 35 3.3.2. Earth-Centered Inertial Coordinate System............................................. 36 3.3.3. Satellite Radial Coordinate System ......................................................... 37 3.3.4. Time ......................................................................................................... 38 3.4. Tracking Systems & Observation Types ........................................................ 39 3.4.1. Deep Space Network (DSN).................................................................... 40 Range .............................................................................................................. 42 Range Rate...................................................................................................... 44 3.4.2. Very Long Baseline Interferometry (VLBI)............................................ 45 3.4.3. Inter-spacecraft Data................................................................................ 48 3.5. Tracking Data Error Sources .......................................................................... 50 x Chapter 4: Orbit Determination .................................................................................. 52 4.1. Orbit Determination History........................................................................... 52 4.2. Classical Orbit Determination......................................................................... 57 4.3. Orbit Determination Problem Formulation..................................................... 59 4.3.1. Aspects of Estimation .............................................................................. 61 4.3.2. Linearization of the Non-linear Problem ................................................ 64 4.3.3. The State Transition Matrix .................................................................... 67 4.3.4. Relating Observations to an Epoch State ................................................ 68 4.4. Orbit Determination Methods......................................................................... 69 4.4.1. Batch Least Squares................................................................................. 69 4.4.2. Weighted Batch Least Squares ................................................................ 72 4.5. Statistical Methods.......................................................................................... 74 4.5.1. Linear Minimum Variance Estimation (LUMV)..................................... 74 4.5.2. Sequential Methods.................................................................................. 76 4.5.3. Comparison of OD Methods.................................................................... 77 4.6. Peturbation Effects.......................................................................................... 79 Chapter 5: Problem Description................................................................................. 81 5.1. Assumptions in the Model .............................................................................. 81 5.2. Batch Processor Algorithm............................................................................. 85 5.2.1. Test Cases ................................................................................................ 85 5.2.2. The Nominal Orbit................................................................................... 86 5.2.3. The Algorithm.......................................................................................... 97 5.2.4. Simulation of Measurement Noise......................................................... 104 5.3. Validity of the Filter...................................................................................... 109 5.4. Orbit Accuracy.............................................................................................. 112 Chapter 6: Results.................................................................................................... 117 6.1. Comments on Orbit Geometry...................................................................... 117 6.2. Measurements Over Time............................................................................. 125 6.3. Main Analysis ............................................................................................... 127 6.3.1. Covariance Analysis .............................................................................. 130 6.3.2. Monte Carlo Runs.................................................................................. 132 6.4. Discussion of Results.................................................................................... 145 Chapter 7: Conclusion.............................................................................................. 150 7.1. Future Work................................................................................................... 152 Appendix A............................................................................................................... 156 Appendix B ............................................................................................................... 157 Bibliography ............................................................................................................. 194 xi List of Tables Table 3.1: Orbital element systems ???.??.????????????.... 34 Table 3.2: DSN DSCC locations [35] ????????????????.?.. 41 Table 5.1: Hughes' nominal orbital elements of LISA & Earth ??????..?.. 89 Table 5.2: Dhurandhar?s nominal orbital elements of LISA & Earth ??..???. 96 Table 5.3: Noise Levels for each measurement [41] ???????????.. 107 Table 6.1: Abridged version of DSN/VLBI/Relative H matrix at day 20 ??.?. 120 Table 6.2: DSN/VLBI H matrix for LISA 1 at day 90 ??????????? 123 Table 6.3: Summary of Mean Absolute Position & Velocity with Standard Deviations ?????...????????. 135 Table 6.4: Summary of Mean Relative Position & Velocity with Standard Deviations ?????????????.... 138 Table 6.5: Summary of Mean Absolute Position & Velocity with Standard Deviations ???????????...??. 142 Table 6.6: Summary of Mean Relative Position & Velocity with Standard Deviations ??????? ????.??. 143 Table 6.7: Parameter Uncertainties (PU) for DSN/VLBI/Relative Case at day 7?149 xii List of Figures Figure 2.1: LISA orbit geometry relative to the Sun and Earth [1] ???????. 7 Figure 2.2: Lasers signals between LISA [1] ???????????????. 11 Figure 2.3: The LISA spacecraft [1] ??????????????????... 12 Figure 2.4: Y-shaped structure [1] ???????????????????.. 12 Figure 2.5: Optical assembly [1] ???????????????????.... 13 Figure 2.6: Orbit geometry of LISA constellation without angles [1] ?????... 18 Figure 2.7: LISA orbit over an entire year [1] ?????????????..?.. 19 Figure 3.1: The two-body problem [47] ???????????????.?? 27 Figure 3.2: Orbital elements [37] ??????????????????.?.. 32 Figure 3.3: Heliocentric-Ecliptic coordinate system (XYZ) [35] ???????.... 35 Figure 3.4: Earth-Centered Inertial coordinate system (ECI) [35] ??????..... 37 Figure 3.5: Satellite Radial coordinate system (RSW) [35] ???..??????.. 38 Figure 3.6: Very Long Baseline Interferometry (VLBI) diagram [42] ????..? 46 Figure 3.7: Pictorial representation of Earth-based DSN and VLBI measurements .. 48 Figure 3.8: Pictorial representation of Earth-based and inter-spacecraft measurements ??????????????????????.. 49 Figure 4.1: The Orbit Determination problem ????????????.??.. 60 Figure 4.2: Along-Track and Cross-Track Errors for LISA [35] ??????.?.. 62 Figure 5.1: Angle measurements to spacecraft ?????????????.?.. 82 Figure 5.2: Geocentric parallax ???????????????????.?. 84 Figure 5.3: Initial position of LISA constellation (Equal Axis) ???????.? 90 Figure 5.4: LISA constellation over 240 days ??????????????? 90 xiii Figure 5.5: LISA position components over 1 year (unperturbed) ??????? 91 Figure 5.6: LISA velocity components over 1 year (unperturbed) ??????? 92 Figure 5.7: Leg length vs. time for an unperturbed orbit over 1 year ?????? 93 Figure 5.8: Percent error in leg lengths for an unperturbed orbit over 1 year ...??. 94 Figure 5.9: Relative velocity for unperturbed orbit over 1 year ???????.? 94 Figure 5.10: Rate of change vs. time of an unperturbed LISA orbit for 1 year ??... 95 Figure 5.11: Dhurandhar magnitude of rate-of-change vs. time for an unperturbed orbit over 1 year ????????????..?.. 97 Figure 5.12: Weighted Least Squares Batch Processor Algorithm ??????? 103 Figure 6.1: LISA constellation at epoch (square axis view) ????????? 118 Figure 6.2: LISA orbit geometry over 20 days ??????????????. 119 Figure 6.3: Eyeball views of triangulation ???????????????.. 122 Figure 6.4: LISA orbit geometry over 90 days ??????????????. 124 Figure 6.5: Earth-based nominal measurements over 1 year ?????????126 Figure 6.6: Inter-spacecraft nominal measurements over 1 year ???????. 127 Figure 6.7: Frequency of observations ?????????????????..129 Figure 6.8: RSS of position sigmas of the Covariance Matrix for LISA 1 ???.. 130 Figure 6.9: RSS of velocity sigmas of the Covariance Matrix for LISA 1 ?.?.? 131 Figure 6.10: MC Simulation, NB, Mean Absolute Position Error ???????. 133 Figure 6.11: MC Simulation, NB, Mean Absolute Position Error with sigmas ??..134 Figure 6.12: MC Simulation, NB, Mean Absolute Velocity Error with sigmas ??.134 Figure 6.13: MC Simulation, NB, Mean Relative Position Error ???????.. 136 Figure 6.14: MC Simulation, NB, Mean Relative Position Error with sigmas ??.. 136 Figure 6.15: MC Simulation, NB, Mean Relative Velocity Error with sigmas ??. 137 xiv Figure 6.16: MC Simulation, noB, Mean Absolute Position Error ??????? 139 Figure 6.17: MC Simulation, noB, Mean Absolute Position Error with sigmas ?? 140 Figure 6.18: MC Simulation, noB, Mean Relative Position Error with sigmas ??..140 Figure 6.19: MC Simulation, noB, Mean Absolute Velocity Error with sigmas ??141 Figure 6.20: MC Simulation, noB, Mean Relative Velocity Error with sigmas ?? 141 Figure 6.21: Frequency of Measurements, Absolute Position Error ??????.. 144 Figure 6.22: Frequency of Measurements, Relative Position Error ?????..?. 145 Figure 6.23: Summary Chart of Days Until Absolute Position Error Equals 1 Kilometer ???????????????????? 146 Figure 6.24: Summary Chart of Days Until Relative Position Error Equals 300 meters ????????????????????.. 147 Figure 6.25: Monte Carlo (Noise, Bias) Parameter Uncertainty of LISA 1 position-components ???????????.????.... 148 xv List of Abbreviations in alphabetical order ALIA Advanced Laser Interferometer Antenna ALIAS Advanced Laser Interferometer Antenna in Stereo ASTROD Astrodynamical Space Test of Relativity Using Optical Devices AU Astronomical Unit BBO Big Bang Observer CHAMP CHAllenging Minisatellite Payload CLIVAR Climate Variability Program cm Centimeter CMNT Colloid Micro-Newton Thruster CSO Compensated Sapphire Oscillator CTE Coefficient of Thermal Expansion DECIGO DECi-hertz Interferometer Gravitational Wave Observatory DOF(s) Degree(s) of Freedom DOR Delta one-way Range DORIS Doppler Orbitography and Radiopositioning Integrated by Satellite DSCC Deep Space Communication Centers CDSCC Canberra Deep Space Communication Center GDSC Goldstone Deep Space Communication Center MS Madrid Deep Space Communication Center DSN Deep Space Network DSS Deep Space Station D/R DSN & Relative data combination xvi DRS Disturbance Reduction System D/V DSN & VLBI data combination D/V/R DSN, VLBI & Relative data combination ECI Earth-Centered Inertial, or Earth-based Geocentric Equatorial (ECI), or IJK Coordinate System EELV Evolved Expendable Launch Vehicle EKF Extended Kalman Filter EM Electromagnetic EOM Equation of Motion ESA European Space Agency FEEP Field Emission Electric Propulsion GEO Geostationary Orbit GFZ GeoForschungsZentrum (i.e. National Research Centre for Geosciences in Pottsdam, Germany) GP-B Gravity-Probe B GPS Global Positioning System GRACE Gravity Recovery and Climate Experiment GRS Gravitational Reference Sensors GTO Geostationary Transfer Orbit Hz Hertz IMS Interferometry Measurement System JPL NASA?s Jet Propulsion Laboratory kg Kilogram xvii km Kilometer lbs Pounds LEO Low Earth Orbit LIGO Laser Interferometer Gravitational-Wave Observatory LISA Laser Interferometer Space Antenna LISAS Laser Interferometer Space Antenna in Stereo LITS Linear Ion Trap Standard LS Least Squares LTP LISA Technology Package LUMV Linear, unbiased, minimum variance m Meter MC Monte Carlo MJ2000 Mean of January 1.5, 2000.0 Reference frame MLE Maximum Likelihood Estimate m/s Meters/second MTPE Mission to Planet Earth NASA National Aeronautics and Space Agency NASA/GSFC NASA Goddard Space Flight Center NASA/JPL NASA Jet Propulsion Lab NB Noise and Bias noB Noise, No Bias NOCC Network Operations Control Center nrad Nanoradians xviii OD Orbit Determination pdf Probability density function P/M Propulsion Module PRARE Precise Range and Range Rate Equipment rad Radians RSW Satellite-Based Radial Coordinate System (where R-axis = Radial, S-axis = Along track, and W-axis = Cross-track) SCI Sun-Centered Inertial, also referred to as Heliocentric (XYZ) sec Seconds SLR Satellite Laser Ranging SMART-2 Small Missions for Advanced Research and Technology STEP Satellite Test of Equivalence Principle STM State Transition Matrix ST7 Space Technology 7 SVD Singular Value Decomposition TCM Trajectory Correction Maneuver TDRSS Tracking & Data Relay Satellite System TRK Tracking DSN Data Type VLBI Very Long Baseline Interferometry WBLS Weighted Batch Least Squares WLS Weighted Least Squares WOCE World Ocean Circulation Experiment xix List of Symbols in order of appearance o Degree s Time-varying strain l Distance ?l Change in distance c Speed of light ?V Change in velocity m 1, m 2 Mass 1 and mass 2 X, Y, Z Position components P Period a Semi-major axis F 1 , F 2 Forces 1 and 2 21 , RR &&&& Accelerations of bodies 1 and 2 21 , RR Distance from origin to masses (m 1 and m 2 respectively) in inertial three- dimensional space )( 21 RRr ?= Relative position vector (i.e. distance between two bodies m 1 and m 2 ) r Magnitude of relative position vector 3 21 3 RRr ?= Magnitude cubed of relative position vector G Universal gravitational constant = 11 3 2 6.6726 0.0005 10 /mkgs ? ? ?? M total Total mass of system (= m 1 +m 2 , ) ? Gravitational constant = GM total xx r r r 3 ?? = && Fundamental equation of motion ? SUN Sun's Gravitational constant (= 1.32712428 x 10 20 m 3 /(solar sec) 2 ) X(t) State vector composed of position and velocity components t Time that represents the State X & ,Y & , Z & Velocity components of the state e r Eccentricity vector (where e is the eccentricity scalar) h r Angular momentum vector (where h is the angular momentum scalar) r r Position vector (where r is the magnitude of position) v r Velocity vector (where v is the magnitude of velocity) r p Periapsis r a Apoapsis i Inclination K ? Unit vector in the K-direction ? Longitude of the ascending node I ? Unit vector in the I-direction n r Line of nodes n i The i th element of the line of nodes vector, n r ? Argument of periapsis e k The k th element of the eccentricity vector, e r ? True anomaly ? Time of periapsis passage ? Mean anomaly xxi E (F) Eccentric Anomaly for an ellipse or circle (comparable quantity for a hyperbola) p Semi-latus rectum ? ~ True longitude of periapsis passage r Orbit radius n Mean motion ? Argument of latitude ? First Point of Aires ? Obliquity of ecliptic (~23.5 o ) I, J, K Three-dimensional axes in Geocentric coordinate system ? Right ascension ? Declination R, S, W Three-dimensional axes on the Satellite Radial coordinate system ? Longitude ? gd Latitude ? Ideal range, magnitude ? v Range vector X sc , Y sc , Z sc Position components of the spacecraft X ts , Y ts , Z ts Position components of the tracking station ZYX uuu ,, Position unit vector components ?& Range rate magnitude xxii ? & r Range rate vector sc X & , sc Y & , sc Z & Velocity components of the spacecraft ts X & , ts Y & , ts Z & Velocity components of the tracking station ij L? Inter-spacecraft range, where i = 1, 2, or 3 and j =1, 2, or 3, and ji ? X i.j , Y i,j , Z i,j Relative position components of the inter-spacecraft legs L 12 , L 23 , L 31 Inter-spacecraft legs ij L?& Inter-spacecraft range rate, where i = 1, 2, or 3 and j =1, 2, or 3, and ji ? ji X , & , ji Y , & , ji Z , & Relative velocity components of the inter-spacecraft legs ? Elevation angle X ? Estimated trajectory (or state) X True trajectory X* Nominal trajectory O ? C Residual Error, or literally "Observed minus Calculated" n-dim n-dimensional column vector of number of state parameters p-dim p-dimensional column vector of observations Y i Column vector of observations (p-dimensional) t i Discrete time )(tX & Time derivative of the state vector (= F(X*,t)) G(X i *, t i ) Observation-state relationship ? i Errors in the observations x State deviation vector, n x 1 y Observation deviation vector, p x 1 xxiii A(t) Partial of the time derivative of the state deviation vector evaluated on the reference trajectory H ~ Observation-state mapping matrix x k Unknown initial state ?(t i ,t o ) State transition matrix I Identity matrix, n x n matrix y i Observation deviation vector f The number of rows in a batch of observations m = p x f, total number of observations H Observation sensitivity matrix J(x) Cost function o x? Least squares correction of the epoch w i Elements of the weighting matrix W Weighting matrix P cov Covariance matrix k x A priori estimate k w A priori weighting (R cov ) i General specified covariance B Unspecified linear mapping matrix, or Baseline in reference to VLBI R Observation noise covariance matrix (where W = R -1 ) K k Kalman gain matrix k P A priori covariance matrix P k A posteriori covariance matrix xxiv ? GEO Geocentric parallax R E Radius of the Earth R E-LISA Distance between Earth and center of LISA constellation L ij Leg length a E Distance from Sun to the Earth Rvs Initial position of center of LISA constellation Vvs Initial velocity of center of LISA constellation M Information matrix N Reduced residual vector i R Position vector i V Velocity components of the State Vector (i.e. [ X & , Y & , Z & ] ) X E , Y E , Z E Position components of Earth M SC Scaled information matrix N SC Scaled reduced residual vector Wsc Scaled weighting matrix SC Scaling Matrix ?T Time-tag error T & ? Clock rate offset ?? Range error due to clock instability, or delta one-way range in reference to VLBI RTLT Round trip light time ? Y (?) Allan Standard deviation ?&? Range rate frequency instability xxv ? Standard deviation, also referred to as sigma RMS Root-Mean-Square RSS Root-Sum-Square, or normal magnitude of the vector norm The normal magnitude of the vector ABSr Absolute OD accuracy in position ABSv Absolute OD accuracy in velocity Li rest Estimated distance to LISA where i =1, 2, 3 Li vest Estimated velocity to LISA where i =1, 2, 3 Li rtru True (or nominal) distance to LISA where i =1, 2, 3 Li vtru True (or nominal) velocity to LISA where i =1, 2, 3 RELr Relative OD accuracy in position RELv Relative OD accuracy in velocity ?L13r, ?L21r, ?L32r Differenced norm of the estimated and true spacecraft positions ?L13v, ?L21v, ?L31v Differenced norm of the estimated and true spacecraft velocities ? ij Correlation coefficients (where i and j equals position (X, Y, Z) and/or velocity (U, V, W) components) U, V, W Also used to indicate three-dimensional velocity components ? 2 Variance ? Number of parameters to be determined from fit PU Parameter Uncertainty 1 CHAPTER 1: Introduction The Laser Interferometer Space Antenna (LISA) mission is a joint effort between the National Aeronautics and Space Administration (NASA) and the European Space Agency (ESA) currently scheduled for launch in 2015 [1]. LISA will function for a five year nominal mission, plus an extended mission of at least 3.5 years once it is placed in its operational heliocentric orbit that trails Earth by about 20 degrees. LISA?s primary objective is to detect and study gravitational waves in the low frequency range of 10 -4 to 10 -1 Hertz (Hz), a region much too low to be perceived by Earth-based observatories. LISA will test Einstein?s theory of general relativity and it is also expected to give further insights into astrophysics. LISA a three-satellite constellation, will measure the relative change in length of adjacent legs of the formation down to picometer (x 10 -12 ) precision. Each leg of the constellation is separated by five million kilometers (km). The leg-length measurement is referenced to a free-floating proof mass that will be shielded from non-gravitational forces inside the optical assembly of each spacecraft using laser interferometry. The relative change in leg length provides information about the direction and polarization amplitude of each gravitational wave. Precise orbit determination (OD) methods are required to initially track LISA and direct each satellite into its proper orbit. During mission execution, accurate OD is required to 2 maintain the appropriate orbit configuration precisely given the extreme sensitivity of scientific data quality to leg length and pointing deviations. OD involves deriving an accurate state vector or initial orbit from which to propagate and map the orbit to a future time and then refining the result to model the real world errors. There are two places in the LISA orbit where OD is particularly important. The first is at the end of cruise phase. At this stage, the Propulsion Modules (P/M) are still actively navigating and correcting the spacecraft orbit thus inter-spacecraft data will not be available. Once the last Trajectory Correction Maneuver (TCM) has occurred and the OD tracking is confident, the propulsion modules can be ejected. The force of the P/M drop-off will give the spacecraft state a "bump." The second crucial OD phase is the commissioning period, roughly a three-month period when LISA will be close to its nominal orbit. Inter-spacecraft lasers will be linked and small adjustments to each satellite?s orbit will be made using micro-thrusters. OD continues to track the spacecraft until the OD is confident and the science of the mission can officially begin. The spacecraft will now be passively controlled using Drag-Free control, and no more large corrections can be made to the orbit during the course of the mission. Hence, OD is crucial to track these spacecraft as they precisely initialize their relative positions, so that the mission can be successful for the scientific operations throughout the nominal and any extended portions. Accurate OD is extremely crucial for successful LISA mission science (i.e. laser interferometry along leg lengths of the constellation) since only passive Drag-Free 3 control will be employed once in operational orbit. The goal of this thesis is to investigate methods to obtain more accurate absolute (with respect to a tracking station on Earth) and relative (with respect to an adjacent spacecraft) OD. This will be done using a weighted batch least squares processor to simulate different combinations and types of observations. The unique aspect of this research is its incorporation of LISA?s inter-spacecraft range and range rate into the OD process that supplements the more standard Deep Space Network (DSN) range, range rate, and directional (right ascension and declination) measurements. Although we study accuracy with the full suite of measurements, we also study accuracy when subsets of all available measurements are used. Such use of "suitable combinations" of observations with a least squares process was studied by researchers as early as Gauss [2]. Currently soft requirements on the mission call for absolute position accuracy to within one kilometer (km) and a relative position accuracy between 300 meters (m) and 20 m. Hechler and Folkner [4] performed a preliminary study using inter-spacecraft data for LISA with its onboard laser tracking system for exclusively the range data type. They found that the incorporation of laser tracking is crucial in determining the arm lengths and achieving tighter accuracy requirements during the experimental phase of the mission, although ground tracking (i.e. range and Doppler) was found to be sufficiently accurate for the transfer and delivery phases of the mission. Hence, this motivated the need to investigate using inter-spacecraft data as a means to improve the OD accuracy. 4 1.1. Thesis Objectives & Contributions This thesis examines the benefits of including relative inter-spacecraft laser data as well as traditional ground tracking data in the OD process to improve both the absolute and relative state estimates. In this thesis, "relative data" is defined as inter-spacecraft range and inter-spacecraft range rate, quantities anticipated to be available from the spacecraft- to-spacecraft interferometry. Earth-based range and Doppler data can be provided by the Deep Space Network (DSN) package, while the angle (i.e. right ascension and declination) data simulates Very Long Baseline Interferometry (VLBI). It is well-known in OD that increasing the amount of data and the time spent observing the target improves accuracy [5]. However, it is not always the case that accuracy improves substantially with the inclusion of new data. This thesis comprehensively evaluates how absolute as well as the relative OD accuracy improves with the inclusion of various measurement combinations. Measurement sets will also be evaluated by their ability to fit within mission constraints and requirements. The effects of longer tracking arcs, greater frequency of measurements, noise, and biases will also be investigated. Our initial hypothesis is that the best OD results will be obtained with the maximum measurement set, including DSN (Earth-based range and range rate), VLBI (angular data), and inter-spacecraft data (inter-spacecraft range and range rate). Specifically, we expect the addition of relative information to yield improvements in both absolute and relative navigation performance. 5 1.2. Overview of Thesis This thesis first overviews the LISA mission in Chapter 2 with a specific focus on LISA?s unique orbit geometry and configuration. During the LISA mission overview, key constraints and technological challenges to this mission will be noted. Related planned and completed missions on which LISA bases elements of its technology will be discussed. Chapter 3 presents the key math and orbital dynamic topics needed in this research. Also, the measurement types, aspects of estimation and tracking data error sources are presented. Chapter 4 describes the orbit determination method used for this work, the weighted batch least squares method. Other methods such as linear, unbiased, minimum variance (LUMV) and sequential methods will be discussed for completeness. A brief OD history is included to provide perspective and further validate choices of methods used in this work. In Chapter 5, the specifics of the nominal LISA orbit, measurement inclusion in the least- squares algorithm, test cases and implementation detail of the simulation are described. The assumptions in the model, validity of the filter, and orbit accuracy test methods will also be defended. Chapter 6 presents the results of the thesis, while Chapter 7 draws conclusions and offers suggestions for future work. 6 Chapter 2: Background The Laser Interferometer Space Antenna (LISA) planned for launch in 2015 is a Cornerstone Mission in ESA?s Cosmic Vision Programme, and is a Great Observatory for NASA?s Structure and Evolution of the Universe 2003 Roadmap: Beyond Einstein Program (Astronomy and Astrophysics Division) [6]. LISA is managed by NASA- Goddard Space Flight Center (NASA/GSFC) and will be operated by NASA-Jet Propulsion Laboratory (JPL). LISA has a five year nominal mission once in operational orbit, plus a possible 3.5 year extended mission [7]. LISA is a constellation consisting of three identical spacecraft, each of which forms a vertex of an equilateral triangle. The distance, or arm length between each spacecraft in the triangle formation, is five million km. The three spacecraft will be in a heliocentric orbit, centered on a reference orbit of one Astronomical Unit (AU), where 1 AU = 149.6 million km from the Sun in the ecliptic plane. The reference orbit center point will trail the Earth by 20 degrees ( o ). Though recently some studies have put this point at 23 o degrees behind Earth [8], the reference point will remain at 20 o in this thesis. There is a tradeoff between choosing 20 o which is more desirable for communications purposes because it is closer to the Earth, while 23 o is better for the drag-free control since it is farther from Earth and its perturbations. The actual orbit will be optimized as the mission draws closer. The entire constellation will be tilted 30 o away from the Sun, at an inclination of 60 o , as illustrated in Figure 2.1. 7 Figure 2.1: LISA orbit geometry relative to the Sun and Earth [1] 2.1. Mission & Science Objectives The scientific objective of LISA is to detect gravitational waves using laser interferometry. In order to measure the changes in distance between two freely-floating proof masses, these masses must be isolated from all non-gravitational forces in each spacecraft [9]. Each spacecraft will track changes in its own motion relative to the other two orbiters [10]. Valuable information concerning the direction of the gravitational wave is gained when comparing the relative motion between two spacecraft of LISA as opposed to following a precise absolute path. LISA aims to be the first space-based detector to identify and study low-frequency (10 -4 to 10 -1 Hz) gravitational waves from three major binary sources: 1) Coalescences of super massive black holes in distant galaxies [1] which are the most violent events in the universe [9] that result from merging galaxies; 2) Stellar-mass black holes spiraling 8 toward massive black holes; and 3) Compact binary star systems in the Milky Way Galaxy [11]. Binaries are simply two objects that orbit each other. Objects can be neutron stars, stars, and black holes for example. Black holes are objects whose gravity is so strong that light cannot even escape it; therefore they cannot be seen, but there is "strong indirect observational evidence" [1] of their existence. Black holes are believed to be formed by the collapse of more massive stars after their nuclear fuel has been expended [1]. This phenomenon emits gravitational wave energy that LISA will be able to detect. A ground-based detector called the Laser Interferometer Gravitational-wave Observatory (LIGO) also views gravitational waves. However, LIGO can only view high-frequency waves from transient phenomena such as supernovae and the final minutes of in-spiraling neutron-star binaries. Low-frequency waves simply cannot be detected with a terrestrial- based detector because of Earth?s changing gravitational field due to atmospheric effects and ground motions [11]. 2.1.1. Gravitational Waves Gravitational waves are vibrations of space and time, or "ripples" in space-time caused by accelerations of masses in space, thought to weakly interact with intervening matter [9]. Gravitational waves are essentially propagating gravitational fields. They create time- varying strain, s, of space-time as defined in Equation 2.1: s = ?l/l (2.1) 9 Where l is a distance and ?l is a change in distance. For LISA, gravitational waves are detected by monitoring the separation of the proof masses. The motion of the source masses can determine further information such as the waveform and the direction to the source. This motion can be inferred by the amplitude and phase variation caused by LISA?s orbital motion [11]. Gravitational waves have been predicted by Einstein?s Theory of General Relativity, but have never been directly detected [9]. In the 1970s, indirect evidence became available to confirm Einstein?s prediction [12]. Einstein?s Theory of General Relativity (1916) proposes matter, space and time are related. The gravity of any mass is the deformation or warping of space-time around it, hence an altering of the path of objects that pass close by. Any variation in movement causes disturbances that become waves in space-time [1]. Mass and gravity are positively correlated: an increase in mass relates to an increase in the gravitational force exerted [13]. LISA aims to test Einstein?s theory of General Relativity and to probe the early Universe. Since gravitational waves are thought to weakly interact with intervening matter, there is a theory that waves generated from the beginning of time would have penetrated the heat and the density of the Big Bang to travel onwards [1]. This could possibly reveal clues of the formation of the early Universe through a search for the continuous spectrum (frequency range) of gravitational radiation generated at the Big Bang [1], the so-called beginning of time for our Universe. The early universe was opaque and not transparent to light during the first 380,000 years following the Big Bang, so light alone cannot 10 reveal the entire history [9]. Thus, the universe is dark and most of what is known about the universe is through electromagnetic (EM) radiation, which interacts with ordinary matter (i.e. gas and stellar material). However, this "dark" matter interacts with gravity; hence gravitational waves give us additional insight [9]. 2.1.2. Laser Interferometry For LISA to detect gravitational waves, laser interferometry precisely monitors the separation of the three spacecraft proof masses. Interferometers are a popular instrument with the Michelson Interferometer named after American physicist Albert Michelson (1852-1931) a common tool used in astrophysics. The famous Michelson-Morley Experiment demonstrated the absence of "ether" a hypothetical medium which would cause the speed of light to be different in different directions. The experiment found that there is no significant motion of Earth relative to ether. This finding led to Einstein?s Theory of Special Relativity in 1906 which states that ether if unobservable simply does not exist. The Michelson-Morley experiment showed that there is no variation in light speed with direction. Laser interferometry takes advantage of the fact that light behaves as a traveling wave that propagates evenly and equally in all directions from its source. Hence, light waves can be identified by their wavelength [1]. Each LISA spacecraft locks the incoming laser signals from the two other LISA spacecraft and transmits them back over the five million km leg length [14]. This is different from the Michelson interferometer, which would reflect the incoming signal. 11 The "phase lock" on the signal ensures that the phase difference is due only to gravitational waves rather than signal loss due to traveling over the substantial leg length and instrument error. It takes about 32 seconds for the laser to transmit a signal to/from spacecraft to spacecraft (~16 light seconds to, ~16 light seconds back) [14]. It is expected that the changes in arm lengths will be measured as accurately as ten picometers [1]. Figure 2.2 illustrates the laser signal paths as they travel between spacecraft. Figure 2.2: Lasers signals between LISA [1] 2.2. LISA Spacecraft Design Each 500 kilogram (kg) LISA spacecraft is a short cylinder (2.7 meter (m) diameter, .89 m height) that supports a Y-shaped tubular structure. The Y-shaped structure encases two instruments, which will be discussed below. Solar panels are mounted on a sunshield that extends from the top of the cylinder. The sunshield forms a cover across the top of the cylinder that will prevent sunlight from hitting the Y-shaped structure and the cylinder walls. Figure 2.3 below illustrates one LISA spacecraft. Note the sunshield and solar panels are not included, so that the inside of the spacecraft can be seen. Also, two 12 X-band 30-centimeter (cm) diameter radio antennae will be mounted on the outside of the spacecraft for communication with Earth. Figure 2.3: The LISA spacecraft [1] The spacecraft and the Y-shaped structure will be graphite-epoxy composite materials because of their low coefficient of thermal expansion (CTE) [11]. The Y-shaped structure is gold-coated both on the inside and outside and is suspended by stressed- fiberglass bands to thermally isolate it from the spacecraft and to reduce radiative heat transfer as shown in Figure 2.4. Figure 2.4: Y-shaped structure [1] 13 2.3. Instrument Package Inside the Y-shaped structure is the optical assembly that contains the optical bench, transmit/receive telescope, and low-power electronics pack [15]. There are two optical assemblies per spacecraft. Figure 2.5 details an optical assembly. Figure 2.5: Optical assembly [1] The optical bench houses the proof mass. Six gold-platinum proof mass cubes (4.6-cm on each side) [10], two per spacecraft, are present in the full LISA constellation. There is a gap between the wall of the structure encasing the proof mass and the proof mass itself. The wall never touches the cube and the spacecraft actively follows the proof mass to maintain the appropriate separation, hence the proof mass continuously "floats" within the spacecraft. For the force on the proof mass to be very low, one needs low stiffness and a low differential displacement of the gap. The two main systems within LISA are: 1) the Disturbance Reduction System (DRS), which includes the Gravity Reference Sensor (GRS) and 2) the Interferometry Measurement System (IMS). The Disturbance Reduction System works to eliminate non-gravitational forces on the proof mass in the measurement direction and to keep the 14 spacecraft pointed in the appropriate direction. The DRS reduces disturbance accelerations to the proof masses to less than 3 x 10 -15 (m/s 2 )/ Hz [3]. Included in the DRS is the GRS, with position sensors to measure the proof mass position, while minimizing coupling to spacecraft motion [16]. Two proof masses are in each GRS. The capacitive sensors also are used to locate the proof mass to ensure that the spacecraft is centered upon the masses. The DRS uses a micro-thrust propulsion system with Micro-Newton (?N) thrusters to perform fine control of spacecraft position and attitude, effectively maintaining a "drag- free" flight condition [16] by commanding the spacecraft to follow the motion of the proof masses. Thrusters must counteract all disturbance forces, such as solar pressure with a position accuracy of 10 -9 m [11]. "Self gravity" between the spacecraft and the proof mass can cause complications in the DRS and create bad measurement noise, the correction of which is an area of current research. In addition, three emerging microthrust propulsion systems are being tested for use on LISA: 1) Field Emission Electric Propulsion (FEEP) using Cesium (Cs-FEEP), or Indium (In-FEEP) propellant, 2) Colloid Micro-Newton Thruster (CMNT), and 3) Precision cold-gas micro-thrusters. The telescope (30-cm. diameter) is used to transmit and receive laser signals from the other spacecraft. The two proof masses are located in each spacecraft 60 o apart from one another. The telescope has an actuator to correct for inter-spacecraft motion, which causes a +/-1 o variation in look-angle to remote spacecraft. The telescope uses interferometer optics of the Interferometry Measurement System (IMS), which measures 15 the change in distance between proof masses to a sensitivity of 4 x 10 -11 m (averaged over one second) [22]. The interferometer measures the changes in the distances between the proof masses separated with leg lengths of five million km in the different spacecraft to picometer-level (x 10 -12 m) accuracy. The inertial sensor, mounted at the center of each optical bench, senses the level of disturbance imposed on each LISA spacecraft. 2.3.1. Drag-Free Control Drag-free control is a concept initially proposed in 1972 and discussed extensively by Lange [17]. To achieve drag-free control, an unanchored proof mass is enclosed in a spacecraft isolating it from external contact forces such as atmospheric drag and solar radiation pressure. Internal forces can be ignored in ideal conditions, thus leaving gravitational forces as the only disturbance. Then the spacecraft can follow the proof mass using low- thrust propulsion. The propagated trajectory is then easier to follow and remain accurate for longer periods of time with only gravity affecting the spacecraft [17], [6]. To be explicit, the external forces from which one must shield the proof mass from include: -Atmospheric drag -Solar radiation pressure -Solar wind -Photon-pressure -Higher order gravitational components -Third-body effects The internal forces also affecting the proof masses including those listed below often cannot be ignored as well: -Self-gravity 16 -Differential gravity -Electric or magnetic field gradients -Radiation pressure 2.4. LISA Launch & Mission Evolution The target launch for the LISA mission is 2015 using an Evolved Expendable Launch Vehicle (EELV) to inject all three spacecraft into an Earth escape trajectory. The spacecraft will separate from the launch vehicle and each other soon after launch then will operate independently. After one orbit, the Earth will be 20 o ahead of the constellation. During cruise phase: a propulsion module (P/M) on each spacecraft makes the appropriate inclination and energy changes in order to rendezvous with final operations orbit [7]. Two circular solar arrays generate the power for the solar-electric propulsion engines. There will be three deterministic maneuvers for each transfer for each spacecraft [7]. However, maneuver execution errors and OD errors will cause the cruise trajectories to deviate from the planned trajectories such that trajectory correction maneuvers (TCMs) must be administered. The propulsion modules will separate from each spacecraft upon delivery to the final orbit (about 13 months after launch). It is expected to require roughly 27 kg (or 59.5 lbs.) of propellant on each spacecraft to achieve the science orbit [1]. The Micro-Newton thrusters that correct for the solar pressure to maintain drag-free control will use a few kilograms of fuel over the mission. 17 Following delivery to operations orbit, each spacecraft will be individually evaluated, calibrated and will begin to operate in drag-free mode during the Element Commissioning Phase. The Acquisition Phase establishes the laser links between spacecraft then the interferometer itself is checked out and calibrated in the Constellation Commissioning phase [7]. The science operation begins after this final commission. The spacecraft positions will evolve under gravitational forces only. The spacecraft will be locally controlled in a "drag-free" manner to keep their positions centered about the proof masses contained within them. 2.5. Orbit Geometry & Configuration LISA is a constellation of three spacecraft in an equilateral triangle configuration with a 60 o (+/- 1 o ) angle between outward-pointing lasers to relative measurements to the other two spacecraft. The LISA constellation is located 20 o behind the Earth. LISA is centered on a reference orbit, a slightly-elliptical Earth-like orbit located one AU from the Sun such that LISA?s period, like Earth?s, is one year. The reference orbit is in the ecliptic plane to minimize transfer change-in-velocity (?V) from Earth [7]. The constellation is inclined 60 o with respect to the ecliptic plane of Earth?s orbit. The 60 o tilt is chosen so that sunlight never enters the interferometer optics; and so that the Sun always illuminates the same part of each spacecraft. Figure 2.6 shows the orbit geometry of the LISA constellation. 18 Figure 2.6: Orbit geometry of LISA constellation without angles [1] Since each spacecraft traces its own path around the Sun, the triangular formation also changes orientation by performing a complete 360 o clockwise rotation once over this same one-year period. Orientation changes are beneficial to help determine the direction of the gravity wave source. The LISA constellation forms an epicycle, or an orbit within an orbit. When viewed from the Sun or Earth, LISA appears to rotate about its center once per year. The rotation is clockwise as viewed from the Sun. Each spacecraft traces an elliptic path as shown in Figure 2.7. 19 Figure 2.7: LISA orbit evolution over an entire year [1] Sweetser [7] describes the orbital motion pattern shown in Figure 2.7, an orbital mechanics phenomenon that only occurs when a plane is tilted 60 o . When a reference point on a circular orbit around a central body is at the center of a plane that is inclined +/-60 o from the orbital plane; all points on the plane maintain a fixed distance from the center to first order. The plane rotates once around its center point per orbit. Dhurandhar et al [18] mathematically validate this theory using the Clohessy-Wiltshire equations. The distance between the spacecraft, referred to as the arm length, is to remain uniform to keep the frequency response of the mission at the chosen optimum [19]. The leg length is five million km and is measured between proof masses. This large separation helps to improve the sensitivity required for detecting low-frequency gravitational waves and requires an appropriately sized laser and telescope. Furthermore, constant distances 20 between the spacecraft minimize problems associated with the inter-spacecraft Doppler shift of the laser signals and allows for optimum frequency response. 2.6. Constraints & Technological Challenges Criteria for choosing the LISA orbits included satisfying at least ten nonlinear constraints imposed by instrument and communication limitations [14]. The constraints are that the leg lengths be kept within +/- 2% of the nominal, the interior angles of the formation must be kept from 58.5 o -61.5 o , the formation should be no more than 20 o behind the Earth for communication and OD requirements, and all three leg lengths must change at a rate less than 15 meters per second (m/s) [14]. LISA's station at 20 o behind Earth is a compromise to: 1) Reduce Earth?s gravitational perturbations on the constellation, such as Earth?s gravity which reduces with distance [7] according to an inverse square relationship, 2) Satisfy launch vehicle constraints, and 3) Enable communication given power requirements as a function of distance from Earth. Two requirements on the operational orbit of LISA are that: 1) the leg length is 5 x 10 6 km +/- 5 x 10 4 km [7], the minimum of which allows for the low frequencies that LISA is designed to detect, and 2) the magnitude of the range rates between spacecraft must meet the 15 m/s constraint described above. These constraints on LISA are present because the leg length range (+/- 5 x 10 4 km) signifies orbit design feasibility. Also, the range 21 rate magnitude specifies angle ranges that onboard interferometers must accommodate [7] throughout the mission (nominal and extended). 2.7. Related Missions Since the 1960?s, attempts have been made to detect high frequency gravitational waves with small interferometers from the ground [20]. The ground-based Laser Interferometer Gravitational-Wave Observatory (LIGO) which was referenced earlier, is a $365-million dollar Earth-based gravitational wave detector located in Louisiana and Washington. LIGO detects gravitational waves at a higher frequency than LISA, including transient phenomena such as supernovae [9]. However, LIGO is currently in need of an upgrade to be Advanced LIGO because detectors were found to experience 100-1000 times more jitter than blueprints indicated. Other ground-based interferometers include Japan?s TAMA, a British and German collaboration GEO600, and an Italy and France partnership known as VIRGO [20]. A set of missions aim to validate subsystems planned for use on LISA and will contribute to the success of future formation-flying interferometer missions. ESA plans to launch the LISA Pathfinder Mission (formerly SMART-2, Small Missions for Advanced Research and Technology, renamed because it will now test technologies that only apply to LISA) in 2009 in a Lagrange point one orbit. LISA Pathfinder will fly two test packages of proof masses and hardware [10]. The payloads will include the complete 22 DRS/GRS payload of NASA?s Space Technology 7 (ST7) to demonstrate the effectiveness of the Gravitational Reference Sensor (GRS) housing two proof masses. This mission will also validate an interferometer design, the Colloid Micro-Newton Thruster (CMNT) by Busek Co., Inc. U.S. [21], and the LISA Technology Package (LTP) [6] also containing two identical proof masses. The LTP drag-free control system consists of an accelerometer (aka inertial sensor), a propulsion system and a control loop [22]. The LTP will instead have the Field Emission Electric Propulsion (FEEP) as thrusters and precision cold-gas micro-thrusters as actuators. The goals of the DRS and LTP are similar, but each uses different technology. LISA-Pathfinder?s mission objectives include demonstrating drag-free and attitude control with two proof masses, assessing the feasibility of laser interferometry at LISA?s accuracy levels, and testing endurance of instrument and hardware in the space environment [22]. The exact design of LISA's proof mass and drag-free control systems will be determined from the studies and feedback of both ST7 and Lisa Pathfinder [22]. Another Beyond Einstein Mission is Constellation-X. Both Constellation-X and LISA are categorized as Great Observatories in NASA?s Beyond Einstein Program, a set of high priority missions that will use X-ray spectroscopy and gravitational waves respectively to study black holes [23]. Constellation-X has a possible launch date in 2011. It will consist of four spacecraft each containing a 1.6 m diameter telescope to measure the spectra of cosmic sources of X-rays [24]. 23 Several other missions have demonstrated technologies applicable to LISA. TRIAD launched Sept. 2, 1972, is the only three-axis drag-free satellite ever flown [25]. TRIAD had a spherical gold-platinum proof mass floating inside of it, unaffected by air resistance and sunlight. Electronic sensors on the spacecraft would fire thrusters if the proof mass approached a wall of the spacecraft to keep it undisturbed from outside forces. The inertial sensors and spacecraft control for LISA are based largely on this mission. It was an experimental spacecraft that flew in a low-altitude orbit until it ran out of fuel. However the spacecraft was tracked for many subsequent years to perform magnetometer experiments. Space Interferometry Mission PlanetQuest (formerly just SIM) developed by NASA-JPL in collaboration with Lockheed Martin Missiles and Space and Northrop Grumman is scheduled for launch in 2011. PlanetQuest will utilize technology related to laser interferometry [26], with a mission of determining the position and distances of stars increased accuracy using optical interferometry. PlanetQuest, like LISA, will also be in an Earth-trailing solar orbit. Gravity Recovery and Climate Experiment (GRACE) consists of two twin spacecraft flying in formation [13]. The precise speed and distance between the two spacecraft are communicated via a microwave K-band ranging instrument. The scientific objective of GRACE is to obtain a highly accurate gravity field map of the Earth. GRACE is located 500 km above Earth. GRACE fulfills goals of NASA?s Mission to Planet Earth (MTPE) - to observe Earth?s gravity field, and responds to many other international organizations 24 concerns such as World Ocean Circulation Experiment (WOCE) and Climate Variability Program (CLIVAR). Gravity-Probe B (GP-B) is a NASA spacecraft launched on April 20, 2004 with a scientific objective to examine two of Einstein?s general theories of relativity: 1) the geodetic effect ? the amount to which the Earth warps local space-time in which it resides, and 2) the frame-dragging effect ? the amount to which the rotating Earth drags local space-time around with it [27]. It orbits the Earth at an altitude of 400 miles. GP-B features nine new technologies, including drag-free control, cryogenics, and new manufacturing and measuring technologies. It is also the first satellite to ever achieve both three-axis attitude control and three-axis drag-free control (the whole spacecraft flies around one of its science gyros while the spacecraft orbits the Earth) [28]. Thus far, GP- B has successfully acquired one-year of science data, but the results are not yet released. CHAMP (CHAllenging Minisatellite Payload) is a German satellite managed by GeoForschungsZentrum (i.e. National Research Centre for Geosciences, GFZ). It is a geopotential research mission that generates highly precise gravity and magnetic field measurements over a five year period. It also uses drag-free spacecraft control [29]. Cassini-Huygens is collaboration between NASA, ESA and Agenzia Spaziale Italiana (i.e. Italian Space Agency) that uses DSN antennas at all three sites continuously [12]. Cassini was launched in 1997, and orbited Saturn on July 1, 2004. Cassini?s probe, Huygens, was dropped onto Titan (a moon of Saturn) in January 2005. Cassini?s 40-day 25 scientific mission began on Nov. 26 2005 and was highly successful. The 40-day mission will be repeated twice over the next two years to vary spacecraft?s position to sensitize the measurements to different gravitational wave directions. Cassini?s overall mission is to measure gravitational waves using radio transmission between Cassini and Earth to reveal velocity changes through the Doppler measurements [30]. Despite not being launched yet, many successors to LISA are also in the planning stages. Astrodynamical Space Test of Relativity Using Optical Devices (ASTROD) by China with arm lengths of 1-2 AU will search for gravitational waves at even lower frequencies than LISA [31]. The National Astronomical Observatory of Japan is proposing the DECi-hertz Interferometer Gravitational Wave Observatory (DECIGO) with shorter armlengths than LISA to bridge the frequency gap between the terrestrial and space- based antennas. NASA and ESA will team together again with Satellite Test of Equivalence Principle (STEP) which will further test "drag free" flight. Another NASA: Beyond Einstein Mission, the Big Bang Observer (BBO) also aims to detect gravitational waves to a finer degree than LISA, perhaps reaching toward data from the beginning of the universe. Lastly, direct extensions to the LISA mission have been proposed with shorter arm lengths and lower noise levels, such as Advanced Laser Interferometer Antenna (ALIA) which will detect intermediate mass black holes [32]. Crowder and Cornish [32] have also proposed the addition of a second constellation in a 20-degree Earth-leading orbit to ALIA (the Advanced Laser Interferometer Antenna in Stereo (ALIAS)) and LISA (Laser Interferometer Space Antenna in Stereo (LISAS)). 26 CHAPTER 3: Methods & Math Techniques This chapter begins with a review of two-body motion, deriving the fundamental equation used to generate an ephemeris from Newton?s Universal Law of Gravitation. Keplerian orbital elements will be overviewed and well as reference frames used in this work. The measurement types and how they will be calculated are introduced. A description of the methods used by Deep Space Network (DSN) to calculate range and range rate is provided along with the method applied by Very Long Baseline Interferometry (VLBI) to calculate angle measurements. The mathematical formulas used to calculate the inter-spacecraft range and range rate will also be shown. The chapter concludes with a discussion of errors present in each the tracking system. 3.1. The Two-Body Problem The OD problem studied in this work was modeled with respect to two-body dynamics since this "simplified" representation is the only gravitational problem for which a closed-form solution has been found. The two-body problem involves two bodies/masses orbiting under mutual attraction in an inertial reference frame, typically with mass m 1 assumed much larger than m 2 such that m 2 does not affect the motion of m 1 . This is illustrated in the three-dimensional (i.e. with dimensions in the X, Y, and Z components) inertial coordinate system in Figure 3.1. 27 Figure 3.1: The two-body problem [47] Note that spherically symmetric bodies, such as planets, behave as point masses. The LISA constellation in this thesis will be treated as a two-body problem between the Sun (the larger mass, m 1 ) and the spacecraft (the negligible mass, m 2 ). The equations of motion for the two masses are derived from Newton?s Universal Law of Gravitation, an inverse square law. Equations 3.1 are based on Kepler?s 3 rd law of planetary motion (i.e. the orbital period is proportional to the semi-major axis to the 3/2 power: P 2 ? a 3 ). () 21 3 21 21 111 RR RR mGm RmF ? ? ? == && (3.1) () 12 3 21 21 222 RR RR mGm RmF ? ? ? == && F 1 and F 2 are the forces and 1 R && and 2 R && are the accelerations of the two bodies, where m 1 and m 2 are the masses of body 1 and 2 respectively. 1 R and 2 R are the distances from the origin to masses m 1 and m 2 respectively in this inertial three-dimensional space, 28 () 21 RRr ?= is the relative position of m 1 with respect to m 2 , while 3 21 3 RRr ?= is the magnitude cubed of the relative position, and G is the universal gravitational constant, which is equal to 11 3 2 6.6726 0.0005 10 /mkgs ? ?? ?. Straightforward manipulation of Equations 3.1 then yields: 21 21 33 3 ()Gm r Gm r Gr m m r rr r ? ??+ =?= && (3.2) Let M total be defined as (m 1 + m 2 ) and ? = GM total , yielding Equation 3.3: .. 3 rr r ? =? (3.3) This gravitational equation of motion (EOM) is dependent only on the masses of the two bodies. The EOM in component form represents a set of three, second-order ( r && ), nonlinear (r 3 term), coupled (X, Y, Z components of r), homogeneous, and autonomous ordinary differential equations [33], where autonomous indicates that a reference epoch for time is not required. The constant ?, the product of G and M total can be determined to a greater accuracy than G or M total alone by accurate Earth satellite tracking [33]. The gravitational constant of the Sun will be defined as, ? sun = 1.32712428 x 10 20 m 3 /(solar sec) 2 . The Equation of Motion will be used to propagate LISA?s trajectory, starting with some initial condition, about its heliocentric orbit for one year as will be seen in Section 5.2.2 (Nominal Orbit) In the least squares problem, the equations of motion must be the exact descriptor of the dynamical system, else dynamical modeling problems may arise [2]. The solution of the EOM requires six independent integrals of constants. The state vector will be defined as 29 an n-dimensional X(t) with dependent variables or constant parameters required to define the time rate of change of the state of the dynamical system [34]. This thesis adopts a typical definition of translational state as a vector of positions and velocities: X(t) = [X Y Z X & Y & Z & ] composed of position components (X, Y, Z) and velocity components ( X & ,Y & , Z & ), and t is the specific time that represents this state. The EOM is used to generate an ephemeris representation of motion using the Cartesian state X(t) or the Keplerian orbital elements. The most accurate way to generate an ephemeris is by using General Relativity to completely describe the EOM by accounting for many perturbations such as effects from the planets and large asteroids in the solar system (mainly the largest ones i.e. Jupiter, Saturn, etc.). A general relativity accuracy level is usually left to describe a planetary ephemeris over a long period of time (i.e. centuries). The numerically integrated EOM based on classical two-body mechanics as described above without the perturbation terms are adequate for most detailed mission analyses and orbit determination and will be presumed sufficient for this thesis. 3.2 Orbital Elements Most of the analysis of this thesis is done using Cartesian elements described above to perform Keplerian transformations defined in Vallado [35], Battin [36] and Chobotov [37]. Orbits are defined by a set of six elements. The first is the eccentricity vector, e r 30 which distinguishes the shape of the orbit and is defined in Equation 3.4. The eccentricity scalar, e, is the magnitude of the eccentricity vector. r rhv e r v r r ? ? = ? (3.4) where h r is the angular momentum vector and is the cross product of the periapsis- pointing position r r and velocity vectors v r , h r =( )vr rv ? , ? is the gravitational constant of the central body. The period, P, of an orbit is derived from Kepler's 2 nd (i.e. equal areas swept in equal intervals of time) and 3 rd (i.e. square of the period of the planet is proportional to the cube of the semi-major axis) laws. The period is defined as: P = ? ? 3 2 a (3.5) where a is the semi-major axis that determines the size, or altitude of an orbit. The semi- major axis can be calculated many ways, from the vis-viva integral or from using the extreme points of the orbit, i.e. the periapsis, r p , and the apoapsis, r a as in Equation 3.6: 2 2 1 2 pa rr v r a ? = ? ? ? ? ? ? ? ? ?= ? ? (3.6) Where, r a = a(1+e) and r p = a(1-e); and r and v are the magnitudes of the inertial position and velocity vectors, respectively. Orbits are described as circular (e = 0), elliptical (e < 1), the hyperbolic (e >1), or parabolic (e =1). Although the LISA orbit is elliptical, some discussion of hyperbolic and parabolic orbits is included for completeness. The inclination, i, defines the tilt of the orbit as seen in Equation 3.7. It is specifically the 31 angle between the unit vector, K ? and the angular momentum vector, h r and varies from 0 o to 180 o . hK hK i r r ? ? )cos( ? = (3.7) The longitude of the ascending node, ? is the angle measured in the equatorial plane positively from the I ? unit vector to the location of the ascending node. The ascending node is where the satellite crosses the equatorial plane from South to North, and varies from 0 o to 360 o . The line of nodes is defined as: hKn r r ?= ? . nI nI r r ? ? )cos( ? =? (3.8) A quadrant check must be employed here that if the i th element, n j , of the line of node vector, n r is greater than 0, then ? = 360 o -?. The argument of periapsis, ?, locates the closest point of the orbit (periapsis) as measured from the ascending node and varies from 0 o and 360 o . en en rr rr ? =)cos(? (3.9) A quadrant check also must be performed on this angle, as it is dependent on the k th - component of the eccentricity vector (i.e., if e k less than 0, then ? = 360 o -?). The last orbital element is the true anomaly, ?, which is the satellite's current position relative to the location of periapsis. re re rr rr ? =)cos(? (3.10) 32 Again a quadrant check is required. If the dot product of the position and velocity vectors is greater than zero, 0 n is satisfied. However, the problem lies in the fact that there are still m unknown observation errors, resulting in m+n total 69 unknowns in only m-equations. The least squares criterion allows a solution to the n- state variables, x, at the epoch t k by providing conditions on the m-observation errors. 4.4. Orbit Determination Methods Orbit determination methods can be used to predict, filter, smooth, and pre-process the observations. Prediction uses existing observations to compute future states via propagation, while filtering is determining the current state using current (and past) observations. Filtering is typically performed by algorithms such as Kalman Filtering and sequential-batch methods that update state and covariance matrices [35]. Smoothing applies to the least-squares and traditional sequential-batch least squares techniques and improves estimates based on existing data to determine the state at a particular epoch. Pre-processing is applying the smoothing process to the observations [35]. Forces acting on each LISA satellite cannot be presumed fully known, but deterministic equations of motion and state transition matrix can be propagated to approximately predict future system states. Measurements are then corrupted by stochastic noise. 4.4.1. Batch Least Squares The least squares criterion as derived by Gauss is applied to approximate the solution of an over-determined (m>n) set of linear algebraic equations. The least squares solution selects a state estimate x that minimizes the sum of the squares of the calculated residuals. 70 As shown in Equation 4.12, x is chosen to minimize the performance index, or cost function J(x). Note ? can be defined from previous Equation 4.11 as ? = y - Hx: ()()HxyHxyxJ T i T i f i T ??==? ? = 2 1 2 1 2 1 )( 1 ???? (4.12) The cost function J(x) is a squared function because a parabola (? 2 ) has a minimum, whereas a line (?) does not [35]. J is also a quadratic function of x, it has a unique minimum when these necessary conditions apply: 1.) ;0= ? ? ? ? ? ? ? ? x J and 2.) 0> ? ? ? ? ? ? ? ? ? ? x x J x x T ?? , for all 0?x? The first condition gives the global minimum/maximum value. The second condition is a second derivative or the Hessian of J(x). The Hessian, which is greater than 0, implies that the symmetric matrix is positive definite indicating that x occurs at a minimum. The first condition can explicitly obtain the normal equation: () ()HxyHHHxy x J T T ??=??== ? ? ? ? ? ? ? ? 0 (4.13) this, can be rearranged to obtain: ( )xHHyH TT ?= (4.14) where x? is the least squares correction that yields the best estimate of x. Equation 4.14 is called the normal equation because if H T H is an m x m symmetric matrix, and if the Hessian is positive definite, then H has full rank. Therefore, the solution for the best estimate of x given the linear observation state relationship Equation 4.11 is Equation 4.15: ( ) yHHHx TT 1 ? ? = (4.15) 71 For the solution in Equation 4.15 to exist, m?n and H must have rank n. If m n (the dimension of y is greater than x) -H is symmetric and positive definite, H has rank = n ()yHHHx TT 1 ? ? = , Least Squares solution for ?x 3.) m < n (more unknowns than equations, but is not typically seen in OD) -H is not square or invertible, H has rank < n ()yHHHx TT 1 ? ? = , minimum norm solution Note that the product H T (HH T ) -1 and (H T H) -1 H T are the pseudo inverses of the H matrix. There are disadvantages to the least squares solution that lead to other methods to solve the OD problem. First, each observation error is weighted equally even though the accuracy of the observations may disagree, an issue addressed with the weighted least squares solution described below. The observation errors may be correlated (i.e. not independent), but the least squares solution does not allow for this. Also, the least squares (LS) method does not consider errors as samples from random processes, i.e. there is no use of statistical information. Due to computational problems that arise from inverting the normal equation matrix (H T H), a variety of solution methods exist for the OD problem. The Cholesky decomposition is a more efficient and sometimes more accurate method to take the inverse, though it is applicable only if the normal equation 72 matrix is symmetric and positive definite [45]. Orthogonal transformation (or the QR factorization method) obtains solutions by applying successive orthogonal transformations to the information array; it yields the same accuracy, but with single- precision computational arithmetic [40]. Given?s rotation, or Householder transformation are two orthonormal transformations that can be used to perform QR factorization [40]. The Singular value decomposition (SVD) method is an alternative to the QR decomposition used to solve the LS problem in a numerically stable manner. SVD is good for detecting and overcoming a possible singularity or near singularity, making it well-suited for ill-conditioned problems. However, the computational effort for the SVD method is more substantial than that for the QR factorization [40]. 4.4.2. Weighted Least Squares The weighted least squares (WLS) method is similar to the least squares except there is an associated diagonal weighting matrix, w i , for each of the observation vectors. The elements of the weighting matrix are normalized with weighting factor values ranging from 0 -1.0. The smaller the weight the less important the measurement (i.e. w i = 0, neglects the observation); while the larger the weighting factor, the more important it is (i.e. w i = 1, gives ?full? weight to that observation). The performance index for the weighted least squares problem is given by: ()() k T kii T i f i T k HxyWHxywWxJ ??==? ? = 2 1 2 1 2 1 )( 1 ???? (4.16) 73 Again, the necessary condition, Equation 4.17 for a minimum J(x k ) requires that the derivative of J with respect to x k be zero: () () k TT k k HxyWHWHHxy x J ??=??== ? ? ? ? ? ? ? ? ? ? 0 (4.17) Equation 4.17 is rearranged to obtain the normal equation: ( ) k TT xWHHWyH = (4.18) Again, if the normal matrix H T WH is positive definite then the solution will have an inverse. The weighted least squares estimate, k x? , is given by: ( ) WyHWHHx TT k 1 ? ? = (4.19) Equation 4.19 can also be expressed as: WyHPx T kk )(? cov = (4.20) where (P cov ) k = (H T WH) -1 . The covariance matrix, (P cov ) k, is a symmetric n x n matrix that must be positive definite to have an inverse. Parameter observability is related to the rank of the (P cov ) k matrix. If all parameters in x k are observable, then (P cov ) k has full rank and an inverse. This coincides with the condition that there is a greater or equal number of independent observations to the number of parameters being estimated (m?n). In general, the larger the magnitude of (P cov ) k , the less accurate the estimate, but caution should be exercised when inferring the accuracy of the estimate by the magnitude of (P cov ) k in the WLS estimate because with normalized weights, (P cov ) k is not valid in the statistical sense [45]. If an a priori value k x is available for x k and a corresponding weighting matrix W k is available, then the following expression can be used for the estimate of x k : 74 ( ) ( ) 1 ? TT kkk k x HWH W HWy W x ? =+ + (4.21) where, k x is the a priori estimate and k W is the weighting matrix for the a priori estimate of k x . 4.5. Statistical Methods The statistical methods are mentioned here for completeness. As seen in the historical overview all of these methods are applications and extensions of Gauss?s least squares method. 4.5.1. Linear Minimum Variance Estimation (LUMV) While the least squares and weighted least squares estimates provide no information on statistical characteristics of measurement errors and no information on a priori errors, the linear minimum variance estimator (LUMV) does address these two issues. Advantages of using the LUMV estimator are that it is relatively simple to use, since a complete statistical description of all the random errors is not necessary and only the first and second moments of the probability density function (pdf) of the observation errors are required (i.e. the mean and covariance matrix associated with random error respectively) [34]. The observation error, ? i , is assumed to be random with zero mean and a specified covariance, (R cov ) i : E[? i ] = 0, E[? i ? i T ] = (R cov ) i (4.22) 75 The state estimation problem is formulated as: x i = ?(t i , t k )x k (4.23) iiii xHy ?+= ~ i = 1, . . ., l where: E[? i ] = 0, E[? i ? j T ] = (R cov ) i ? ij and ? ij = Kronecker delta function. The best linear, unbiased, minimum variance estimate, k x? , of the state x k is shown below. The linear condition implies that the estimate is made up of a linear combination of the observations: Byx k =? (4.24) where, B is an unspecified linear mapping (n x m) matrix selected to obtain the best estimate [34]. The unbiased condition follows: E[ k x? ] = x = E[By] = E[BHx k + B? k ] = x k Recall that the expected value, or mean of the error is zero: E[? k ] =0, so that: BHx k = x k , to arrive at Eqn 4.25: BH = I (4.25) Lastly, the minimum variance condition is the task of selecting the (n x m) matrix B to minimize the covariance matrix. The covariance matrix, (P cov ) k for an unbiased estimate is: (P cov ) k = E[( k x? - E[ k x? ])( k x? - E[ k x? ]) T ] With substitutions of the linear and unbiased conditions and unspecified Lagrange multipliers, L, the covariance matrix, (P cov ) k is given by: (P cov ) k = (H T R -1 H) -1 . And the LUMV estimate of x k is Equation 4.26: k x? = (H T R -1 H) -1 H T R -1 y (4.26) 76 Note that if the weighting matrix, W, is equal to the inverse of the observation noise covariance matrix, R, i.e. W= R -1 , the LUMV solution agrees with the WLS solution of Equation 4.21. Similarly to WLS, the LUMV estimate can be modified to propagate the estimate of the state x k and its associate covariance matrix. The LUMV estimate with a priori information can be determined from these propagated values [34]. The LUMV estimate of Equation 4.26 also will agree with the maximum likelihood estimate (MLE) developed by Fisher (1912) if the observation errors are assumed to be normally distributed with a zero mean and covariance, (R cov ) k [47]. 4.5.2. Sequential Methods Gauss implied that observations can be inaccurate and of unknown and unknowable errors [2], leading to modern statistical methods such as the Kalman filter. Sequential estimators process observations as soon as they are received. The solution and its associated covariance then can be mapped to another time [47]. A major advantage of this method is that it avoids some of the matrix inversion problems of the least squares methods, in that the matrix to be inverted will be the same dimension as the observation vector [45]. The estimate of the state x k , is simply Equation 4.27, which is obtained by recursively solving the supporting equations [ ] kkkkkk xHyKxx ?+=? (4.27) where, k x = ?(t k ,t j ) j x? K k = weighting matrix (aka Kalman gain matrix) = [ ] 1? + k T kkk T kk RHPHHP P k = a priori covariance matrix = ?(t k ,t j )P j ? ? (t k ,t j ) 77 P k = a posteriori covariance matrix = ( T k HR k -1 H k + P k -1 ) -1 A disadvantage of the sequential processor is that the state estimate covariance matrix will approach zero (P k ? 0) as the number of observations becomes large. The gain K k also approaches zero, and the estimation procedure becomes insensitive to observations and will diverge due to errors that have now been introduced into the linearization procedure, computational errors, and errors due to an incomplete math model [45]. To improve on this disadvantage, two modifications can be used: 1) the Extended Sequential estimator (aka Extended Kalman-Bucy Filter), which minimizes the effects of errors due to the neglect of higher order terms in the linearization process; and 2) the State noise compensation method to combat errors in the dynamic model. Additional modifications such as the 1) Q-Matrix compensated sequential estimation algorithm; and 2) The Dynamic Model Compensation (DMC) estimation algorithm, a sequential estimation method that compensates for the unordered effects in the differential equations which describe the motion of an orbiting vehicle [34]. 4.5.3. Comparison of OD Methods The batch least squares method is often compared to the sequential methods to justify their use, as these are the two basic data processing procedures for tracking data used to estimate a spacecraft's orbit. A disadvantage of both the batch and sequential processors is that if the true and reference state are not close together, the linearization assumptions leading to Equation 4.6 are not valid and the estimation process will diverge [45]. The 78 batch and sequential algorithms are both iterated until convergence, while the Extended Kalman Filter (EKF) will obtain near convergence in a single iteration [45]. Batch processing allows for simultaneous use of multiple measurements, hence making it most useful for post-flight analysis of very large observation datasets. It provides an estimate of the state at some chosen epoch using an entire batch or set of data. The estimate and its associated covariance are mapped to other times. Process noise in the batch processor complicates the solution of the normal equations, the dimension of the normal matrix increases from n (#of state parameters) to m (# of observations) [47]. Without process noise, the sequential algorithm can be shown to be mathematically equivalent to the batch process. Given the same data set, both will produce the same estimates when mapped to the same times [47]. Process noise is added in the EKF however, and this algorithm is not equivalent to the batch processor, but has been shown to be very close [47]. The numerical integrator in the EKF is restarted at each observation time (good for real-time applications), and some state noise is needed to ensure that convergence does not occur. State noise compensates for various error sources in the processing of ground-based or onboard data [47]. As noted before, the covariance matrix, (P cov ) k , may approach zero as the number of observations becomes large, causing filter divergence, though process noise is usually added to deal with this issue. Sequential Filters are best suited for real-time estimation when the mathematical model used to describe the system is accurate, though corrections for unordered force effects can also best be implemented with sequential estimation algorithms [34]. 79 4.6. Perturbation Effects Major sources of error can be divided into two main groups: satellite force model effects and measurement model inaccuracies. Error sources are mainly satellite-dependent, meaning that the satellite configuration, orbit altitude, inclination, and measurement systems used for each individual mission plays a major role in the effects of errors. Errors affecting the satellite force model are further subdivided into gravitational and nongravitational parameters. Gravitational parameters include: mass of the Earth (GM), Geopotential coefficients (C lm , S lm ), solid Earth and ocean tide perturbations, mass and position of the moon and planets, and general relativistic effects. Observations of range/range rate and angular data from various satellites work to improve Earth?s gravitational model [47] so that the gravitational potential of Earth, a major error source for LEO satellites, can be better compensated. With deep space orbits like that used for LISA, each satellite is mainly affected by solar-radiation pressure and third-body effects. Gravitational forces from large planets (i.e. Jupiter, Saturn) and asteroids may pose a greater error source on LISA than the gravity of Earth. Non-gravitational parameters include drag, solar and Earth radiation pressure, thrust and other effects such as magnetic forces. This research presumes a two-body force model in which each satellite is only influences by the Sun's gravitational effects. Appendix A summarizes basic calculations that verify third body effects may be presumed negligible for this preliminary OD work. Errors in the measurement model include: inertial and terrestrial coordinate system errors and ground-based measurements errors. Inertial and terrestrial coordinate system errors are precession, nutation, and polar motion. Ground-based measurement errors include 80 tracking station coordinates, atmospheric effects, instrument modeling, tectonic plate motion and clock accuracy [47]. Other appreciable error sources include inaccuracies in the mathematical model of the dynamical (satellite force) and instrument measurement system. Computer errors are due to storage and execution time problems or errors in data processing and the fact that various parameters, such as ?, in the math model are only approximately known are all additional sources of error. 81 CHAPTER 5: Problem Description This chapter describes the LISA OD simulation that was developed for this thesis. Following a summary of assumptions, the batch processor algorithm will be presented, including the four data type cases tested, the initial condition of the LISA orbit, and data noise sources. Methods used to identify errors in the algorithm and software and a discussion of the equations required to obtain the absolute and relative OD errors are presented. 5.1. Assumptions in the Model The major modeling assumptions in the simulation are: 1) The values for parameters ? SUN and any other constants are considered to be known exactly. 2) Earth and the LISA spacecraft are all considered to be point masses. 3) The tracking station(s) are located at the center of this point mass Earth. The rotation of the Earth and DSN tracking stations are not taken into account. A specific Deep Space Station antenna is not modeled. 3a.) DSN range and range rate measurements will be two-way tracking made from the single tracking station on the point mass Earth 82 3b.) VLBI measurements will also be made from that single tracking station, despite the practical requirement for two stations (Section 3.4.2). Ultimately, the angle measurements are the output as illustrated in Figure 5.1. To avoid confusion VLBI refers to the angular measurements right ascension, ?, and declination, ?. Figure 5.1: Angle measurements to spacecraft 4) The system is unperturbed and will be treated as a two-body system with the Sun being the central body force and each spacecraft in the LISA constellation the negligible mass (see Appendix A for a summary of error resulting from this assumption). 4a.) The placement of Earth in the model is considered known and is 20 degrees ahead of the reference center of the LISA constellation. The Earth is also a negligible mass in the two-body system. ? Y Z X ? ? 83 4b.) The reference orbit is deterministic, a system whose time evolution can be predicted exactly. Therefore, Newtonian acceleration will be integrated and the trajectory of the spacecraft and Earth will be propagated with the Equations of Motion (Section 3.1) for one year from some arbitrary epoch t = 0. 5.) Tracking will be based on X-band uplink/downlink Doppler and range to data to the spacecraft. 5a.) Noise and bias for DSN and VLBI will be based on values for the year 2000 X-band system for DSN using typical error values as listed in Cathy Thornton's monograph [41]. 5b.) Values for inter-spacecraft range and range rate random error and biases were made to match equivalent Earth-based range and range rate values. 6) The random noises added to the actual measurements will be assumed to follow a Gaussian/Normal distribution. 7) A three month (or ~90 day) commissioning period at the end of cruise phase, when LISA is close to its nominal orbit, is the key timespan that this thesis will investigate. During this crucial phase, P/M drop-off has occurred and OD tracking continues until drag-free control and the science of the LISA mission can officially begin. Calculation of the geocentric parallax, ? GEO , roughly justifies assumption #3, or the angular difference between observations from the center of Earth and those from a site [35]. A simple, not to scale, drawing as seen in Figure 5.2 illustrates geocentric parallax. 84 Figure 5.2: Geocentric parallax The geocentric parallax, ? GEO , is: ? GEO = arctangent(R E / R E-LISA ) = arctangent( 6,378,135 meters / 5.1955 x10 10 meters) = 1.2276 x 10 -4 radians, or ~0.007 degrees The radius of Earth is 6,378,135 meters, negligible compared to the distance from Earth to the reference center of LISA which is 5.1955 x 10 10 m. Hence, the approximation of the tracking station being located at Earth's center is reasonable, since the geocentric parallax is less than a tenth of a degree (= 0.007 degrees). For assumption # 4, recall that the initial set-up of LISA placed 20 o behind the Earth was a compromise to be sufficiently far to reduce the gravitational pull of the Earth and Moon on this constellation, yet still be ?close? enough to meet radio communication requirements with Earth, and to meet launch vehicle propellant requirements. While drag free control shields the proof mass from perturbations such as solar flux, the gravity of other planets still does affect the spacecraft. However, back of the envelope calculations of the force effects of the two major perturbation sources, Jupiter at closest approach to Sun, and Earth at closest approach to LISA. Appendix A uses Newton's Law ? GEO R E R E-LISA 85 of Gravitation to compare these effects to the force of the Sun on LISA to justify the assumption of a two-body system. It was found that the Sun poses the most force on LISA in both cases, which supports the fact that this can be considered a two-body problem for now. Certainly a full perturbation analysis is required prior to launch. 5.2. Batch Processor Algorithm The computational algorithm for the batch processor used in this thesis is based on the Batch Processing algorithm from Tapley's text Statistical Orbit Determination (2004) [47]. The MATLAB 7.0 code that implements the entire simulation and batch processing algorithm is presented in Appendix B. The LISA OD problem requires modeling three spacecraft and four main test cases described below. 5.2.1. Test Cases Four test cases that utilize different combinations of observation data were investigated: Case 1: DSN & VLBI (D/V) Case 2: DSN, VLBI, & Relative (D/V/R) Case 3: DSN only (DSN) Case 4: DSN & Relative (D/R) Note that cases such as Angles-only, Relative-only and Angles/Relative observation cases were also considered, however, they were eliminated due to the fact that system state is not fully observable with theses measurement sets. Case 3 explores accuracy with strictly the DSN data package, which includes range and range rate measurements. Case 86 1 uses a standard DSN and VLBI (angles) package. Case 2 incorporates the more novel relative data (i.e. inter-spacecraft range and range rate) available from the LISA spacecraft into DSN and VLBI and is expected perform the best according to the hypothesis that more data types would equal better OD. Case 4 combines DSN and relative data measurements. 5.2.2. The Nominal Orbit Once a nominal orbit is determined for the baseline and extended missions, proper OD is needed to achieve and maintain the LISA constellation in the specified configuration. An approximate knowledge of the orbit is required for least squares analysis in order to supply the initial "guess" of the state to the algorithm. Section 4.3 stated that the true state of the spacecraft in the OD problem is never known. However, these results were simulated without actual data. It is assumed in this simulation that the truth orbit is known and is the nominal orbit; which will be stated below. The test of how well the weighted batch least squares (WBLS) algorithm estimates is based on how well the final estimate agrees with the truth. Noises and biases will be added into the measurements to simulate real world data obtained from instruments and tracking stations. Nominal orbits for use on the actual LISA mission have been determined through the works of Steve Hughes [14] and S.V. Dhurandhar et al [18]. Hughes? nominal orbit is found to maintain a nearly equilateral triangle configuration for two-body orbital dynamics and to be unstable when all the planets and Earth's moon are added to the 87 system. It is also a good guess for numerical optimization as an appropriate orbit can be determined for LISA in a perturbed orbit. In [14], Hughes creates five performance measures, or cost functions for his orbit optimization problem. The challenge is to "find solutions that minimize a cost function and also satisfy a set of non-linear constraints imposed by instruments and communication limitations." There are ten nonlinear solution constraints for LISA: 1) Leg lengths must be +/-2% of the nominal 5 x10 9 m (3 constraints for each of the 3 legs) 2) Interior angles of the formation should remain between 58.5 o -61.5 o (3 constraints) 3) The formation must orbit no more that 20 o behind the Earth to meet communication and OD constraints (1 constraint) 4) The maximum average rate-of-change must be <15m/s for any of the legs (3 constraints) Considering the solution constraints, Hughes numerically optimizes the performance of the formation in the presence of perturbations for these five cost functions (cf n , where n= 1, 2, . . ., 5) formulated based on science and/or practical engineering or OD concerns: cf 1 : Difference in rate-of-change between two legs in the formation (science). cf 2 : Difference in lengths of two legs of the formation (science). cf 3 : Average rate-of-change of all three legs of the formation (engineering). cf 4 : Maintaining equal leg lengths over time (science). 88 cf 5 : Deviation of interior angle(s) of the formation from 60 o as an alternative method to ensure that all three leg lengths are equal (science). Hughes improves the nominal orbit, by using it as an initial guess through his optimization approach for each of the five cost functions and checks their feasibility by how well they satisfy the constraints. Since this thesis assumes an unperturbed orbit, Hughes? orbit will be used, though Dhurandhar?s will be commented upon. In the nominal orbit, four orbital elements (a, e, i, ?) are the same for all orbits in the formation. The longitude of the ascending node, ?, and mean anomaly, ?, are different for each orbit. Below is a summary of the nominal orbital elements from [14]: Nominal orbit: a = a E = 1 AU = 149597870691 m e = E a L 32 = 9.6484 x 10 -3 , where L = leg length = 5 x 10 9 m i = E a L 2 = 1.6711 x 10 -2 radians = 0.9574956 o ? = ?/2 or 3?/2 (90 o , or 270 o ) for LISA 1: (?, ?) for LISA 2: (? + 2?/3, ? ? 2?/3), where 2?/3 = 120 o for LISA 3: (? - 2?/3, ? + 2?/3) P = 1 year = ? ? 3 2 r = 3.1558 x 10 7 sec. ~ 365.26 days where ? Sun = 1.327124 x 10 20 m 3 /s 2 The first four elements (a, e, i, ?) are the same for all spacecraft. The right ascension, ?, and mean anomaly, ?, have a phasing of 120 o . The periods, P, of each LISA are all one year (or 3.1558 x 10 7 sec) - the same as the Earth's period. Table 5.1 below 89 summarizes the thesis initial condition for all three LISAs and Earth at an epoch unattached to any time system. Table 5.1: Hughes' nominal orbital elements of LISA & Earth at epoch = 00:00:00 a ( m ) e i ( o ) ? ( o ) ? ( o ) ? ( o ) LISA 1 1.4960 x 10 11 0.0096484 0.95750 0 90 0 LISA 2 1.4960 x 10 11 0.0096484 0.95750 120 90 -120 LISA 3 1.4960 x 10 11 0.0096484 0.95750 240 90 120 Earth 1.4960 x 10 11 1.3438 x 10 -16 0 NaN NaN 0 The right ascension, ?, and the argument of periapsis, ?, of Earth are unobservable, indicated as "NaN" (i.e. Not a Number) because the inclination is 0. The work of this thesis was done using Cartesian elements. The nominal orbit derived in Table 5.1 was transformed to a position and velocity using Orb2X2.m (see Appendix B). The Cartesian state can be transformed back to the Keplerian orbital elements as a check using X2Orb.m (see Appendix B). The position and velocity of the reference center of the LISA orbit was derived to be: Rvs = [ 0 a E 0] Vvs = [-v 0 0] Where a E is the semi-major axis that is initially all in the +Y-direction (Heliocentric); and v is the tangential velocity of the reference orbit and is instantaneously in the negative X- direction. The Earth was rotated to be 20 o ahead of the reference center, as illustrated in an "equal-axis" view of Figure 5.3, the initial position for the nominal orbit of the three LISA spacecraft, Earth and the Sun (in RSW coordinates). The initial position and velocities of Earth and LISA are then used as the initial condition to propagate the orbit using the LISA_EOMderivs.m (see Appendix B). 90 -6 -4 -2 0 2 x 10 10 0 0.5 1 1.5 2 x 10 11 -5 0 5 x 10 9 W : C r o ss- T r ack, m SH: Initial Condition Positions of LISA S: Along-Track, m R: Radial, m Sun LISA 1 LISA 2 LISA 3 Earth Reference Center Figure 5.3: Initial position of LISA constellation (equal axis view) The motion of the constellation around the Sun is carried out for 240 days using the unperturbed equations of motion as illustrated in Figure 5.4. Figure 5.4: LISA constellation over 240 days Note that each LISA rotates around its epicycle and is phased 120 o from the other spacecraft as was indicated in the initial orbital elements of Table 5.1. This is clearly -2 0 2 -2 -1 0 1 2 x 10 11 -5 5 x 10 9 S: Along-Track, m LISA trajectory over 240 days R: Radial, m W: Cross-Track, m Sun LISA 1 LISA 2 LISA 3 Earth 91 shown in Figures 5.5 and 5.6 which illustrate the sinusoidal motion of the position and velocity components of each LISA. 0 1 2 3 4 x 10 7 1.4045 1.4046 1.4046 1.4047 1.4047 1.4047 x 10 13 Time (sec) N o r m o f Po s i ti o n (m ) Norm of Position vs. Time 0 1 2 3 4 x 10 7 -2 -1 0 1 2 x 10 11 Time (sec) X- c o m pone nt of P o s i tion ( m ) X-Component vs. Time 0 1 2 3 4 x 10 7 -2 -1 0 1 2 x 10 11 Time (sec) Y - c o m pon e n t o f P o s i t i on ( m ) Y-Component vs. Time 0 1 2 3 4 x 10 7 -3 -2 -1 0 1 2 3 x 10 9 Time (sec) Z- c o m pon e n t o f P o s i t i on ( m ) Z-Component vs. Time LISA 1 LISA 2 LISA 3 Earth Figure 5.5: LISA position components over 1 year (unperturbed) 92 0 1 2 3 4 x 10 7 2.7963 2.7963 2.7964 2.7965 2.7965 2.7965 x 10 6 Time (sec) N o r m of V e lo c i t y ( m /s e c ) Norm of Velocity vs. Time 0 1 2 3 4 x 10 7 -4 -2 0 2 4 x 10 4 Time (sec) X - co mp o n en t o f V e l o ci t y ( m / s e c ) X-Component vs. Time 0 1 2 3 4 x 10 7 -4 -2 0 2 4 x 10 4 Time (sec) Y - c o m p one nt o f V e loc i t y ( m /s e c ) Y-Component vs. Time 0 1 2 3 4 x 10 7 -500 0 500 Time (sec) Z- c o m p one nt o f V e loc i t y ( m /s e c ) Z-Component vs. Time LISA 1 LISA 2 LISA 3 Earth Figure 5.6: LISA velocity components over 1 year (unperturbed) The leg length, L ij , is a scalar defined as Equation 5.1: jiij rrL ?= where i?j and i & j = 1, 2, 3 (5.1) The initial leg lengths for the nominal orbit are: [L 12 L 23 L 31 ] = [5.0093 x 10 9 4.9752 x 10 9 5.0093 x 10 9 ] meters Figure 5.7 illustrates the sinusoidal leg length variation over an unperturbed orbit propagated for one year. 93 0 0.5 1 1.5 2 2.5 3 3.5 x 10 7 4.975 4.98 4.985 4.99 4.995 5 5.005 5.01 5.015 5.02 5.025 x 10 9 Time (sec) Le g Le ngt hs ( m ) SH: Leg Lengths vs. Time Leg 12 Leg 23 Leg 31 Figure 5.7: Leg length vs. time for an unperturbed orbit over 1 year Figure 5.8 is the differenced leg length graph in which the nominal leg length of exactly 5 x10 9 m is subtracted from each leg length. The result is a graph that shows the percent error. The leg length varies by +/- 0.5 %, or +/- 25,000 km. The velocity differences or relative velocities between the legs are shown in Figure 5.9. 94 0 0.5 1 1.5 2 2.5 3 3.5 x 10 7 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 x 10 7 Time (sec) Le g Le ngt hs D i f f e r e nc e ( m ) SH: Leg Lengths Difference vs. Time Leg 12 Leg 23 Leg 31 Figure 5.8: Percent error in leg lengths for an unperturbed orbit over 1 year 0 0.5 1 1.5 2 2.5 3 3.5 x 10 7 4 5 6 7 8 9 10 x 10 -3 Time (sec) M a g n it u d e o f R e lat i ve V e lo ci t y ( m / s ) SH: Relative Velocity vs. Time Leg 12 Leg 23 Leg 31 Figure 5.9: Relative velocity for unperturbed orbit over 1 year 95 The magnitude of the rate of change of the leg length is given in Equation 5.2, which is a derivative with respect to time of the position vector. dt rd dt rd dt Ld j i ij ?= (5.2) The rate-of-change of the LISA orbit vs. time is depicted in Figure 5.10. Note that the mission constraint is satisfied when the rate-of-change does not exceed 15 m/s. 0 0.5 1 1.5 2 2.5 3 3.5 x 10 7 -4 -3 -2 -1 0 1 2 3 4 Time (sec) M a gni t ude of R a t e - o f - C h a nge ( m / s ) SH: Magnitude of Rate-of-Change vs. Time Leg 12 Leg 23 Leg 31 Figure 5.10: Rate of change vs. time of an unperturbed LISA orbit for 1 year The Hughes' nominal orbit is calculated in NomOrb.m in Appendix B. Though ultimately not used in this thesis, the nominal orbit of Dhurandhar [18] is shown below. The successful operation of LISA requires the stability of the LISA constellation and constant leg lengths between spacecraft in a perturbed system. Dhurandhar et al. explicitly demonstrate how any configuration of the spacecraft lying in the planes making +/- 60 degree angles with the ecliptic using the Clohessy-Wiltshire equations remains stable, in 96 that the inter-spacecraft velocities remain constant to first order with respect to the dimensions of the configuration. The trajectory using Dhurandhar's equation is shown in Table 5.2. Table 5.2: Dhurandhar?s nominal orbital elements of LISA & Earth a ( m ) e i ( o ) ? ( o ) ? ( o ) ? ( o ) LISA 1 1.4963 x 10 11 0.009681 0.96673 0 90 0 LISA 2 1.4969 x 10 11 0.0094387 0.95631 118.92 88.349 -117.27 LISA 3 1.4969 x 10 11 0.0094387 0.95631 241.08 91.651 117.27 Earth 1.4960 x 10 11 1.3438 x 10 -16 0 NaN NaN 0 Dhurandhar proposed using the Clohessy-Wiltshire equations to find a proper orbit for LISA. His initial orbit does not stray too far from nominal guess orbit that Hughes proposed. The argument of periapsis, ?, is around 90 o for all three spacecraft and the phasing on ? and ? are roughly 120 o . The inclinations are small (<1 o ) and the eccentricity is nearly zero, for a nearly perfectly circular orbit. There was some variation found for the period P of the orbit using Dhurandhar's equations: P (Earth) = 3.1558 x 10 7 sec P (LISA 1) = 3.1567 x 10 7 sec P (LISA 2 & 3) = 3.1587 x 10 7 sec Figure 5.11 reveals why the Dhurandhar orbit was not used in this thesis. For the unperturbed orbit, Dhurandhar's orbit does not remain constant. Figure 5.11 illustrates the magnitude of rate-of-change extending beyond the 15 m/s mission constraint. In future work, Dhurandhar's orbit would be used in the perturbed system. 97 0 0.5 1 1.5 2 2.5 3 3.5 x 10 7 -100 -50 0 50 100 150 Time (sec) M a gni t ude of R a t e - o f - C h a nge ( m / s ) DHUR: Magnitude of Rate-of-Change vs. Time Leg 12 Leg 23 Leg 31 Figure 5.11: Magnitude of rate-of-change vs. time for Dhurandhar's unperturbed orbit over 1 year 5.2.3. The Algorithm The goal of the batch processor is to estimate the state deviation vector x o at a reference time t o in order to give the best estimate of the state at that time. Time t o is an arbitrary epoch, and all quantities are assumed to be mapped to the epoch using state transition matrix ?(t i ,t o ). The input values include initial conditions of reference state X*(t o ) and optional a priori state estimates. Although available in the algorithm, the use of a priori estimates will not 98 be addressed here. Equation 5.3 is the normal equation form of o x? for the weighted least squares. Recall that W=R -1 : ( ) WyHxWHH T o T =? (5.3) Equation 5.3 is iterated until o x? , the least squares correction, converges. Note the solution of o x? requires the inversion of information matrix M: ( )WHHM T = (5.4) The reduced residual vector, N, is the right-side of Equation 5.3: WyN T ?= The M and N matrices are accumulated over the batch of measurements to keep the dimensions of the matrix from growing and to allow for easier inversion as shown below in Equations 5.5 and 5.6: 1. [] ),( ~ ),( ~ 1 oiii T f i oii T ttHWttHWHH ??= ? = (5.5) 2. [] ii T f i oii T yWttHWyH ? = ?= 1 ),( ~ (5.6) The state transition matrix ?(t i ,t o ) is attained from integration of ()()( ) kk tttAtt ,, ?=? & . A(t) is the partial of the time derivative of the state vector F(X * ,t), evaluated on the reference trajectory X*: ( ) X tXF tA ? ? = ? , )( (5.7) The code that calculates the state transition matrix (LISA_STMderiv.m) is listed in Appendix B. The observation-state mapping matrix i H ~ is shown below: 99 X tXG H ii i ? ? = ? ),(~ (5.8) where G(X i * , t i ) is the observation-state relationship evaluated on the nominal, or reference trajectory. The partials are computed by taking the partial derivatives of the observation with respect to the reference trajectory. A listing of the partials is below and will give insight into which measurements give position and velocity information: Range Partial Derivatives Range with respect to position: ? ? ? ? ? ? ? ? ??? = ? ? i Ei i Ei i Ei i i zzyyxx R ??? ? )()()( Range with respect to velocity: ? ? ? ? ? ? ? ? = ? ? 000 i i V ? Range Rate Partial Derivatives Range Rate with respect to position: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? i i i Ei i i i i Ei i i i i Ei i i i ZZ Z YY Y XX X R ? ? ? ? ? ? ? ? ? ? & & & & & & & Range Rate with respect to velocity: ? ? ? ? ? ? ? ? ??? = ? ? i Ei i Ei i Ei i i zzyyxx V ??? ? )()()(& Right Ascension Angle Partial Derivatives Right Ascension with respect to position: ( ) ( ) () ? ? ? ? ? ? ? ? ? ? ? = ? ? 0 coscossin 2 iEiE ii i i XXXXR ???? Right Ascension with respect to velocity: ? ? ? ? ? ? ? ? = ? ? 000 i i V ? 100 Declination Angle Partial Derivatives Declination with respect to position: () () ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? )cos( )sin(1 tantan i i i i i i i i i i i i i i Z YX R ? ? ? ? ? ? ? ? ? ? ? Declination with respect to velocity: ? ? ? ? ? ? ? ? = ? ? 000 i i V ? Inter-spacecraft Range Partial Derivatives Inter-spacecraft Range with respect to position (LISA i range from LISA j): ()()( ) ? ? ? ? ? ? ? ? ??? = ? ? Lij ji Lij ji Lij ji i Lij ZZYYXX R ??? ? Inter-spacecraft Range with respect to velocity: ? ? ? ? ? ? ? ? = ? ? 000 i i Lij V ? Inter-spacecraft Range Rate Partial Derivatives Inter-spacecraft Range Rate with respect to position: () ( ) () ( ) () ( ) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? +? ? +? = ? ? Lij Lij ji Lji Lij Lij ji Lji Lij Lij ji Lji i Lij ZZ ZZ YY YY XX XX R ? ? ? ? ? ? ? ? ? ? 131313 & && & && & && & Inter-spacecraft Range Rate with respect to velocity: ()()( ) ? ? ? ? ? ? ? ? ??? = ? ? Lij ji Lij ji Lij ji i Lij ZZYYXX V i ??? ?& As shown above, the inter-spacecraft range will be computed with respect to the other spacecraft rather than to the Earth tracking station (denoted X E , Y E , Z E ). For the angle 101 partial derivatives, MATLAB was used to compute and verify the equations in a compact form. Appendix B contains the MATLAB code (H_partials.m) that calculates the partial derivatives. The observation sensitivity matrix, H i , is given by: ),( ~ oiii ttHH ?= (5.9) When inverting the accumulated M, information matrix and the N reduced residual matrix, a question of observability arises. Observability is crucial for the matrix inversion to work (i.e., sufficient data must exist for all states to be known), since the covariance matrix (P cov ) is needed to solve the normal equation (i.e. (P cov ) = M -1 ). The inverse of a matrix can only be taken if the matrix is square (m x m) and nonsingular. However, not all square matrices have an inverse. In order to be nonsingular, one of these conditions musts be satisfied: 1) no row (column) is a linear combination of the other rows (columns), or 2) the determinant of the matrix is not equal to zero [49]. The first condition is related to the rank of a matrix. Only a matrix composed of independent equations will be full rank, where rank(M) = m. Therefore, the rank check is a quick way to determine observability. Each navigation measurement combination case had to be treated separately in order for the entire program to run cohesively. The M and the N matrices were sized depending on the case. For the cases without relative measurements, the initial M is a 6x6 matrix of zeros, and the N is a 6x1 vector of zeros, and a separate M and N were accumulated for 102 each spacecraft. For the cases with relative measurements, the additional data requires M to be an 18x18 vector while N is an 18x1 vector. For the relative cases, all three spacecraft states were solved simultaneously. However, in order to obtain observability in all cases, scaling matrices had to be added so that the inverse of M could be computed. The accumulation of M and N is described as: M SC = M SC + ((H i *SC) T *W*(H i *SC)) N SC = N SC + ((H i *SC) T *W*(y i )) where, SC is the scaling vector. The scaling vector, SC is a square diagonal matrix that does not scale position terms but multiplies velocity terms by 1x10 -6 . SC has the same dimensions as M. In order to remove the scaling vector numerical effects in the results, another scaling matrix was multiplied into the normal equation: SCSCSCSC NPSCNMSCx *)(***? cov 1 == ? or: NPNMx **? cov 1 == ? X o * is held constant for the beginning of each iteration. The initial guess orbit, X o * is augmented by the least squares correction value, o x? , after each iteration: nonono xXX )?()()( 1 += ? ?? (5.10) The recursive nature of the algorithm continues until convergence occurs, when the WBLS algorithm outputs the best estimate of the state, X ? , which is equivalent to the final corrected X o * . Presented in Figure 5.12 is the weighted batch least squares batch processor algorithm used to estimate the LISA orbit, it is based on an algorithm in [47]. The batch processor code (LISA_WBLS.m) and the code that ran the entire process (MCRUNbp.m) is listed in Appendix B. 103 Figure 5.12: Weighted Least Squares Batch Processor Algorithm B1. Integrate Reference Trajectory & State Transition Matrix from t i-1 ? t i - ( )( )ttXFX , ?? = & , initial condition X*(t i-1 ) ?X*(t i ) - ? ? ? ? ? ? ? ? ? = X tXF tA ),( )( -()( ) ( ) oo tttAtt ,, ?=? & , initial condition ?(t i-1 , t o ) ? ?(t i , t o ) B. Read in Time, Observation (Obs.), Weights t i , Y i , W i B2. Accumulate Current Observations - ? ? ? ? ? ? ? ? ? = X tXG H i i ),(~ -y i = Y i ? G(X i * , t i ) - ),( ~ oiii ttHH ?= - ii T i HWHMM += and ii T i yWHNN += Yes No Else Update & Iterate: X o * ? X o * + o x? XSTOP ? ? Has Process Converged ? tolx o 30) however, for a normally distributed sample the larger value of n, the better the approximation [50]. This allowed for statistics such as the mean and standard deviations of the outputs to be obtained. Monte Carlo is a very time consuming approach because of the number of cases needed to obtain 114 meaningful statistics, but it does have the capture the effects of non-linearities in the systems; whereas, the covariance analysis is a less time-consuming process, but it ignores any non-linearities in the system. The summary graphs of Chapter 6 will plot the absolute and relative OD error of position and velocity. The absolute OD accuracy is divided into the position and velocity absolute accuracy, represented by the root-sum-square (RSS), or norm expression: RSS = () )()( 2 21 2 21 2 21 zzyyxxnorm ?+?+?= The absolute OD accuracies in position (ABSr) and velocity (ABSv) are taken as the norm of the RSS difference between the estimated LISA state (i.e. L1r est and L1v est ) and the nominal LISA state (i.e. L1r tru and L1v tru ) as shown below. The relative position (RELr) and velocity (RELv) error of the constellation is obtained by first taking the norm of the spacecraft estimated leg length (L13, L21, L32) minus the norm of the corresponding truth data leg length: then the norm of all three legs are calculated as shown below. [ ]( ) []())33(),22(),11( )33(),22(),11( truesttruesttruest truesttruesttruest vLvLnormvLvLnormvLvLnormnormABSv rLrLnormrLrLnormrLrLnormnormABSr ???= ???= 115 []() []vLvLvLnormRELv rLrLrLnormRELr vLvLnormvLvLnormvL vLvLnormvLvLnormvL vLvLnormvLvLnormvL rLrLnormrLrLnormrL rLrLnormrLrLnormrL rLrLnormrLrLnormrL trutruestest trutruestest trutruestest trutruestest trutruestest trutruestest 3221,13 3221,13 )23()23(32 )12()12(21 )31()31(13 )23()23(32 )12()12(21 )31()31(13 ???= ???= ???=? ???=? ???=? ???=? ???=? ???=? The current inter-spacecraft accuracy requirement of 300 meters [8] may be tightened in the future. Anticipating this constraint, a stricter relative position error of 20 meters will be used in this thesis. The mean position and velocity error of each LISA is also outputted, along with the final batch processor estimated state vector ( X ? ). This value is used to compute and compare leg length differences and orbital elements. The sample variance,? 2 is expressed in Equation 5.17: i T i Wyy m ? ?? = 1 2 ? (5.17) Where, i T i Wyy ? = sum over all observations (recall that y i is the observation deviation vector and W is the weighting matrix) m = total number of Measurements = p x f ? = number of parameters to be determined from the fit The parameter uncertainty, PU, for each state element can be obtained from the sample variance as shown in Equation 5.18: )( cov 2 PdiagPU ?= ? (5.18) 116 where, the diag(P cov ) term refers to the diagonal elements of the covariance matrix, P cov corresponding to the state, which was also outputted. The parameter uncertainties are added or subtracted from the estimated state to give a final answer with uncertainties. Other important outputs include RMS of each LISA and number of iterations. The period and semi-major axis are two key orbital elements that will also be assessed. They are related through: P = 2?* ? 3 a . There are mission constraints on both elements. Maximum period error is 38 seconds from the nominal [8]. Period accuracy is especially crucial since timing affects the range rate measurement. The semi-major axis should be within the range of 0.9 AU and 1.1 AU [8], which indicates that 0.1 AU is the maximum semi-major axis error since 1AU is the nominal. The weighted batch least squares algorithm is expressed in the LISA_WBLS.m (see Appendix B) script which calls the aforementioned scripts for the nominal spacecraft and Earth tracking station trajectory propagated using the EOM, the STM, the simulated observations, and the partials of the observations. MCRUNbp.m (see Appendix B) runs the entire process for all test cases and generates comparative statistics in a Monte Carlo analysis. 117 CHAPTER 6: Results This chapter will provide an overview of the orbit geometry of this problem. Next, it will consider the observation sensitivity matrix (H) and the measurements, which will give insights into the main results that follow. All simulation results and some figures were generated in MATLAB. 6.1. Comments on Orbit Geometry A closer look at the LISA simulation orbit geometry is illustrated in the following Figures at the epoch time (= 0), 20 days, and 90 days, respectively. Twenty days is where DSN- only becomes effective, assuming a good initial guess. Ninety days is roughly the end of the assumed three-month commissioning period required to put LISA in the proper orbit, so that the drag-free control and the scientific portion of the mission can commence. Figure 6.1 is a square axis view of the set-up at epoch. 118 Figure 6.1: LISA constellation at epoch (square axis view) Over LISA?s one year orbital period, the three spacecraft rotate into positions of zero Doppler (i.e. range rate measurement) when viewed from Earth. At these points, the range and Doppler only reveal information about the degrees of freedom (DOFs) along the line to the Earth; hence, the problem is not observable. This explains why DSN-only (i.e. range and Doppler range rate) measurements are not able to track LISA until a significant arc has been observed. Mathematically, observability is when the M matrix is full rank. This is found to occur around day 20 when constellation has sufficiently rotated to enable more visibility of all three spacecraft -6 -5 -4 -3 -2 -1 0 1 x 10 10 0 0.5 1 1.5 2 x 10 11 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 x 10 9 Z (m) Initial Condition (SH) Positions of LISA X (m) Y (m) Sun LISA 1 LISA 2 LISA 3 Earth Reference Center 119 from Earth as seen in Figure 6.2 below (the shape indicates position of the body at day = 20 and the line is the path traveled). Figure 6.2: LISA orbit geometry over 20 days The addition of Relative data (i.e. inter-spacecraft range and range rate) adds information in the directions not observed directly by range and Doppler. For the relative cases, the state of the entire constellation is estimated at once from inverting the 18x18 information matrix, M. Extra information is obtained about the legs of the constellation, and the location (i.e. range) and velocity (i.e. range rate) of the spacecraft with respect to each other. -8 -6 -4 -2 0 x 10 10 0 5 10 15 x 10 10 -1.5 -1 -0.5 0 0.5 1 1.5 2 x 10 9 X (m) LISA trajectory over 20 days Y (m) Z (m) Sun LISA 1 LISA 2 LISA 3 Earth 120 The M matrix is an accumulation over the measurement batch and contains the observation sensitivity, or H matrix. The partials in the observation-state mapping matrix H ~ give observability information on each measurement and multiplying it by the STM maps the observability to a specific time. An abridged version (18x6) of the H partial matrix in Table 6.1, for the DSN/VLBI/Relative case at day 20 takes an explicit look at how the measurements affect the determination of the estimated state, for LISA 1. Table 6.1: Abridged version of DSN/VLBI/Relative H matrix at day 20 DSN/VLBI/Rel. LISA 1 LISA 2 LISA 3 Day =20 Position Velocity measurement X Y Z X & Y & Z & ?. ?. ? L1 8.8E-01 4.7E-01 4.6E-02 0 0 0 0 0 1L?& -9.6E-08 1.8E-07 -2.9E-09 8.8E-01 4.7E-01 4.6E-02 0 0 ? L1 -7.9E-13 -4.3E-13 2.0E-11 0 0 0 0 0 ? L1 -9.4E-12 1.7E-11 0 0 0 0 0 0 ? L13 3.3E-01 -4.0E-01 8.5E-01 0 0 0 0 X 13L?& -9.8E-08 -1.7E-08 3.0E-08 3.3E-01 -4.0E-01 8.5E-01 0 X ? L2 0 0 0 0 0 0 X 0 2L?& 0 0 0 0 0 0 X 0 ? L2 0 0 0 0 0 0 X 0 ? L2 0 0 0 0 0 0 X 0 ? L21 -6.1E-01 -5.6E-01 5.6E-01 0 0 0 X 0 21L?& -3.6E-08 -9.3E-08 -1.3E-07 -6.1E-01 -5.6E-01 5.6E-01 X 0 ? L3 0 0 0 0 0 0 0 X 3L?& 0 0 0 0 0 0 0 X ? L3 0 0 0 0 0 0 0 X ? L3 0 0 0 0 0 0 0 X ? L32 0 0 0 0 0 0 X X 32L?& 0 0 0 0 0 0 X X In Table 6.1, each measurement contributes to knowledge of the state at a different level. Range (e.g., ? L1) from Earth only provides information most prominent in the X- and Y- 121 components, whereas, the Z-component is less discernible because the LISA constellation as viewed from Earth is smallest along the Z-axis. Similarly, the range rate from Earth (e.g., 1L?& ) has the largest partials (x 10 -1 ) in the X & - and Y & -component of velocity, and it reveals much smaller information about the Z & - component. Range rate also gives some position information, albeit, of much lesser value than for velocity. However, for the inter-spacecraft data, range (? L13 and ? L21) discloses information of equal magnitude (x 10 -1 ) of all three position components. Moreover, inter-spacecraft range rate (e.g., 13L?& and 21L?& ) provides data for all three velocity components. Declination (e.g., ? L1) is most observable in the Z-component (x 10 -11 ), a magnitude which appears small compared to range, or range rate, but is large for angles measured in the units of radians. Right ascension (e.g., ? L1) is most observable in the X and Y-component. Table 6.1 clearly shows that LISA 1 is sensitive to leg 13, but also gives information about leg 21. The same trend is seen for LISA 2, which is set up to measure range and range rate of leg 21, but also contains information on leg 32. LISA 3 includes information on leg 32 in addition to leg 31, where the "X" in Table 6.1 denotes a value for that measurement in the abridged chart. The 60 degree angle of the constellation triangle means that any one spacecraft is well-observed in the plane of the constellation by the other two; since information can still be obtained about the location of the spacecraft despite poor DSN data through this triangulation. Figure 6.3 illustrates this concept of triangulation and how, with the addition of relative range and range rate data, 122 more information is obtained about a spacecraft. In the illustrations, an "eyeball" is ranging to one LISA (in the order LISA 1, LISA 2 LISA 3), but picks up relative information about the adjacent two legs. Figure 6.3: Eyeball views of triangulation Figure 6.3 also supports an earlier statement (Chapter 5) that the measurement noise reduces as the square root of the number of measurements. A complete relative inter- spacecraft dataset is similar to having a single spacecraft dataset measured three times as LISA 2 LISA 1 LISA 3 LISA 3 LISA2 LISA 1 LISA 2 LISA 1 LISA 3 123 often. Thus, relative data should effectively be the same as lowering the noise of the range and Doppler data types by a factor of 3. In the non-relative cases, the states of the three spacecraft are estimated separately and only information pertaining to that spacecraft is measured, as seen in the H matrix (6x6) below for LISA 1 for DSN/VLBI at day 90. Table 6.2: DSN/VLBI H matrix for LISA 1 at day 90 DSN/VLBI LISA 1 Day = 90 X Y Z X & Y & Z & ? L1 -1.61E-01 9.87E-01 1.69E-04 0 0 0 1L?& -2.02E-07 -3.29E-08 -1.01E-08 -1.61E-01 9.87E-01 1.69E-04 ? L1 5.51E-16 -3.39E-15 2.04E-11 0 0 0 ? L1 -2.01E-11 -3.27E-12 0 0 0 0 For both H matrices (Tables 6.1 and 6.2), there is no velocity partial information attained from the range (?), inter-spacecraft range (?Lij), or angle measurements (? and ?). Velocity information is achieved only with range rate (?& ) and inter-spacecraft range rate ( Lij?& ). All measurement partials contain position information. However, there is no Z- component information in right ascension angle (?). The components of Table 6.2 retain the same signatures as discussed in Table 6.1, though the direction and size are dissimilar at day 90 due to the different orbit station. Angle data is shown to help initially by augmenting unobservable DSN, which will be evidenced by obtainable results of DSN/VLBI for small tracking arcs. With a good initial guess, DSN/VLBI values can be achieved for a tracking arc of one day. The angle data given with VLBI is the right ascension (?), which locates the spacecraft in the X-Y 124 plane, and declination (?), which orients the spacecraft vertically. Knowledge of distances, or positional components are provided by this data. Long arcs of Doppler and range can provide angular position accuracy of 40 nrad or better [41]; hence this could explain why the angular data with accuracy of 50 nrad is not useful compared to range and Doppler after a sufficiently long tracking arc has been traversed. Figure 6.4 shows the LISA simulation at the end of the three-month commissioning period. Due to the significant tracking arc all spacecraft are observable from Earth. Figure 6.4: LISA orbit geometry over 90 days -15 -10 -5 0 x 10 10 0 5 10 15 x 10 10 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x 10 9 X (m) LISA trajectory over 90 days Y (m) Z (m) Sun LISA 1 LISA 2 LISA 3 Earth 125 6.2. Measurements Over Time Figures 6.5 and 6.6 examine how measurements from Earth, and the inter-spacecraft measurements, would look over the course of a year. Note the low declination of the LISA orbit, varying from +/- three degrees. The range and range rate follow a sinusoidal pattern. The inter-spacecraft range and range rate varies in a sinusoidal fashion, hence, recall these values were plotted earlier in the nominal orbit discussion (Section 5.2.2). The inter-spacecraft range varies about its nominal leg length (i.e. 5 x 10 9 m). A mission requirement is that the LISA arm length should be +/- 1% of the 5 x 10 9 m [8] (i.e. between 4.95 x 10 9 m and 5.05 x 10 9 m), even with noisy and biased data, inter-spacecraft range values fall within that constraint in this two-body problem set-up. 126 0 100 200 300 400 4.9 5 5.1 5.2 5.3 5.4 5.5 x 10 10 Time (days) R a nge ( m ) LISA Range vs. Time 0 100 200 300 400 -600 -400 -200 0 200 400 600 Time (days) Ra n g e Ra t e ( m / s ) LISA Range Rate vs. Time 0 100 200 300 400 -200 -100 0 100 200 Time (days) R i ght A s c e n s io n A n gle ( d e g r e e s ) LISA Right Ascension Angle vs. Time 0 100 200 300 400 -3 -2 -1 0 1 2 3 Time (days) D e c l i n at i o n A n g l e ( d eg r e es ) LISA Declination Angle vs. Time LISA 1 LISA 2 LISA 3 Figure 6.5: Earth-based nominal measurements over 1 year 127 0 50 100 150 200 250 300 350 400 4.97 4.98 4.99 5 5.01 5.02 5.03 x 10 9 Time (days) I n te r - S/ C R a n g e (m ) LISA Interspacecraft Range vs. Time LISA 1 LISA 2 LISA 3 0 50 100 150 200 250 300 350 400 -4 -2 0 2 4 Time (days) I n te r - S/ C R a n g e R a te (m / s ) LISA Interspacecraft Range Rate vs. Time Figure 6.6: Inter-spacecraft nominal measurements over 1 year 6.3. Main Analysis The main analysis was done on the four main cases: 1) DSN/VLBI (D/V) 2) DSN/VLBI/Relative (D/V/R) 3) DSN only 4) DSN/Relative (D/R) Note that the "Relative" refers to the inter-spacecraft range and range rate data. For the main analysis, a series of 35 randomly-generated simulations for each case were executed 128 to perform a Monte Carlo analysis. The means and standard deviations of outputs were collected to represent error statistics. The cases were run for the roughly three-month (or 90-day) commissioning period. A summary of the Monte Carlo analysis characteristics are listed below, given noises and biases previously described in Section 5.2.4: - Random number generated noise - Bias (Yes and No) - Random generated Initial Guess Error if indicated - Each case was run from 1 to 90 days - Ran 35 times per day per case - With 1 measurement/day Two Monte Carlo runs were performed on each case: 1) Noise, Bias, Random Initial Guess Errors, 1 measurement/day (NB) 2) Noise, No Bias, No Initial Guess Error, 1 measurement/day (noB) One case of three measurements per day was also done on the noise, bias Monte Carlo (NB) run to compare the effects of more frequent measurements on the estimation. Prior to this main analysis, single runs were also performed over the entire orbital period (i.e. 365 days) to observe trends. Since each case was run just one time, the single runs were faster than the Monte Carlo runs. A random number generator was not used in the single runs; so in some cases, the noise was turned on at its full value, which is treated as an additional bias since it is now constant; and other cases turned off the noise. While single runs are not the focus of the analysis, they did allow data patterns to be observed, over longer time periods, with larger biases, as well as testing more measurements per day. The covariance run was done on a case with no noise, no bias, and no initial guess error to generate a covariance plot on each case. In essence, the covariance run reveals the 129 standard deviation of each case. Recall from Chapter 5 that the covariance matrix is a square and symmetrical matrix with the variance (or standard deviations squared) along its diagonal, which discloses the quality of the state estimate. The covariance run described below presumed the case: 1) No Noise, No Bias, 1 measurement/day with - No Initial Guess Error - Each case was run from 1 to 365 days - Each case ran 1 time per day Figure 6.7 illustrates the frequency of observations obtained at each update interval over the entire tracking arc. Figure 6.7: Frequency of observations Obs. Site Obs. Site Obs. Site Update Interval 24 hrs. (8 hrs.) Observations (Obs.) 1 meas./day (3 meas./day) Tracking Arc (Fit Span) days (90 or 365 total) Satellite Pass over Observation Site Time 130 6.3.1. Covariance Analysis A no noise, no bias run is really a test of how well the filter works. This case appropriately computes 0 for all the error outputs. The diagonals of the covariance matrix and the estimated state are the only information extracted from this case. The first insight seen early in testing is that the DSN-only case was found to be unobservable until day = 20, therefore, its tracking starts then. For this set-up, it took two iterations for all cases to converge to their tolerance. The root-sum-square (RSS) of the square root of the covariance matrix (P cov ) diagonals were plotted, i.e. the standard deviations (sigma, or ?). Figure 6.8 shows the RSS of position sigmas for LISA 1. The covariance graphs for LISA 2 and 3 show similar trends. Position Sigma RSS of Covariance Matrix for LISA 1: No Noise, No Bias Single Run 0 1 10 100 1000 10000 100000 1000000 10000000 1 7 20 35 45 60 75 90 11 0 14 0 17 0 20 0 23 0 26 0 29 0 32 0 35 0 Tracking Arc (days) P o s i t i on S i gm a R S S ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.8: RSS of position sigmas of the Covariance Matrix for LISA 1 Around day 30, DSN/VLBI/Relative and DSN/Relative, overlap and it takes about 90 days for DSN/VLBI and DSN to overlie. These overlaps suggest that VLBI or angle data 131 eventually becomes ineffectual compared to the other measurement types. DSN/VLBI/Relative starts out with the least standard deviation, followed by DSN/VLBI. This suggests the importance of VLBI observations early in the mission phase. At 15 days, DSN/Relative surpasses DSN/VLBI. After day 15, the two cases with Relative information have the least standard deviation for the remainder of the run, concluding that Relative data gives a clear advantage in accuracy over non-relative data. Figure 6.9 illustrates the RSS of velocity sigmas for LISA 1: Velocity Sigma RSS of Covariance Matrix for LISA 1: No Noise, No Bias Single Run 0.00000001 0.0000001 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 1 7 20 35 45 60 75 9 0 1 10 1 40 1 70 2 00 2 30 2 60 2 90 3 20 3 5 0 Tracking Arc (days) V e l o c i t y Si g m a R SS ( m / s ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.9: RSS of velocity sigmas of the Covariance Matrix for LISA 1 In Figure 6.9, the DSN/Relative velocity sigma surpasses DSN/VLBI much earlier by day 3, as indicated in the graph. DSN/VLBI and DSN overlap at day 70. DSN/VLBI/Relative and DSN/Relative overlie around day 30. In both Figures 6.8 and 6.9, there is a decreasing trend to all lines. As the tracking arc and number of measurements increase, the RSS sigmas of the covariance decrease. 132 DSN/VLBI and DSN mimic each other; and DSN/VLBI/Relative overlaps DSN/Relative. These figures indicate that with Relative data VLBI adds value for the first 30 days. Without Relative data, it is advantageous to have VLBI for a longer time as indicated in both graphs judging by how long the DSN/VLBI case outperforms DSN. The addition of Relative data improves the overall performance, as well. While VLBI helps initially, it is the presence of Relative data that gives a clear advantage in spacecraft state knowledge for the remainder of the orbital period. 6.3.2. Monte Carlo Runs The Monte Carlo (MC) runs were performed for three months (or 90 days), the commissioning period as well as a period considered sufficient based on previous runs to capture most measurement/convergence dynamics. Each case was iterated 35 times per day. The noises were randomly generated. The initial guess error was only randomly generated for the Monte Carlo (noise/bias denoted "NB") case. From day 1 to 19, the initial guess error was randomly generated with a maximum of +/-10 km in position and +/-1 x 10 -5 km/sec in velocity. However, after 20 days when the DSN-only case was added the random generator was reduced to 1 km since DSN would often be "close to singular, or badly scaled," thus, unable to converge with even that much initial guess error. No checks were done on the upper limit of how much guess error the other cases could take, yet they were consistently able to converge with the prior stated offset (+/- 10 km) offset. In addition to DSN requiring a sufficiently long tracking arc, the initial guess for DSN cannot be too far from the actual state. For any case, a bad initial guess could 133 lead to divergence of the solution. However the other cases can compensate better for a bad initial guess because they have additional observation types (angles and/or relative). Figure 6.10 shows the mean absolute position error. For this simulation, random guess errors were random +/- 10 km, therefore, the DSN measurement was only observable from approximately day 30. Thus, its tracking began there. Monte Carlo: Mean Absolute Position Error: Noise & Bias, Random Initial Guess Error Case 1 10 100 1000 10000 100000 1000000 10000000 100000000 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M ean A b sol u t e P o si t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.10: MC Simulation, NB, Mean Absolute Position Error The 1-sigma standard deviation error bars of Figure 6.11 show that DSN/VLBI/Relative remains distinct until day 15 then it begins to overlap with DSN/Relative. Around day 75, DSN/VLBI/Relative overlaps with DSN/VLBI. By day 10, DSN/Relative and DSN/VLBI are overlapping. DSN remains distinct from all the other data indicating it will never be as accurate over a 90 day tracking arc. Figure 6.12 shows the mean absolute velocity error with 1-sigma (1-?) standard deviation error bars. 134 Monte Carlo: Mean Absolute Position Error with 1-sigma error bars: Random Noise & Bias, Random Initial Guess Error Case 1 10 100 1000 10000 100000 1000000 10000000 100000000 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M ean A b s o l u t e P o si t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.11: MC Simulation, NB, Mean Absolute Position Error with sigmas Monte Carlo: Mean Absolute Velocity Error w/ 1-sigma error bars: Noise & Bias, Random Initial Guess Error Case 0.00001 0.0001 0.001 0.01 0.1 1 10 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M e an A b sol u t e V e l o ci t y E r r o r ( m / s ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.12: MC Simulation, NB, Mean Absolute Velocity Error with sigmas Table 6.3 is a summary of the raw data of mean absolute position and velocity errors and 1-sigma standard deviations for the cases at 1, 20, 30, 60, and 90 days. There is the expected 3 (~ 1.7) improvement for absolute OD accuracy of the Relative cases. The 135 absolute error improvements are much greater than 3 for shorter tracking arcs, though the gap between the Relative and non-Relative cases decrease over time. Table 6.3: Summary of Mean Absolute Position & Velocity with Standard Deviations Mean Absolute Position Error (m) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 5294.6 (+/- 1639.9 ) 2293 (+/- 1195) ----- 4.26 x 10 6 (+/- 3.00 x 10 6 ) 20 1743.2 (+/- 699.57) 300.99 (+/- 224.83 ) ----- 881.07 (+/- 665.41) 30 1426.6 (+/- 493.73) 273.03 (+/- 150.03) 1.60 x 10 7 (+/- 7.15 x 10 6 ) 363.43 (+/- 169.2) 60 632.86 (+/- 229.53) 254.07 (+/- 20.436 ) 2.30e+05 (+/- 97931) 255.49 (+/- 20.851) 90 188.92 (+/- 86.714) 112.54 (+/- 4.7928) 19378 (+/- 11069) 112.61 (+/- 4.786) Mean Absolute Velocity Error (m/s) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 0.066 (+/- 0.030) 0.00257 (+/- 0.00193 ) ----- 0.16246 (+/- 0.10881) 20 0.0017 (+/- 6.7 x 10 -4 ) 3.19 x 10 -5 (+/- 2.16 x 10 -5 ) ----- 7.24 x 10 -5 (+/- 5.11x 10 -5 ) 30 8.93 x 10 -4 (+/- 3.35 x 10 -4 ) 4.04 x 10 -5 (+/- 1.52 x 10 -5 ) 2.635 (+/- 1.273 ) 4.88 x 10 -5 (+/- 1.64 x 10 -5 ) 60 1.5 x 10 -4 (+/- 9.30 x 10 -5 ) 5.79 x 10 -5 (+/- 3.92 x 10 -6 ) 0.049 (+/- 0.025 ) 5.81 x 10 -5 (+/- 3.97 x 10 -6 ) 90 4.69 x 10 -5 (+/- 2.10 x 10 -5 ) 4.43 x 10 -5 (+/- 1.29 x 10 -5 ) 0.0049 (+/- 0.0018) 4.43 x 10 -5 (+/- 1.29 x 10 -6 ) Note that on day 90, there is still a standard deviation of 11 km for the DSN absolute position error, which is a lot of uncertainty in the answer. While the Relative cases have a mean absolute position of 112 m, which is not too far from 188.92 m for DSN/VLBI, the standard deviation on the relative cases is far less (= +/- 4.79 m). Thus far, the 136 Absolute data is shown to improve with the presence of Relative data. Figure 6.13 shows that the relative error also improves with the presence of Relative data. Monte Carlo: Mean Relative Position Error: Random Noise & Bias, Random Initial Guess Error Case 1 10 100 1000 10000 100000 1000000 10000000 100000000 1 3 5 7 101520253035374045505560657075808590 Tracking Arc (days) M ean R e l a t i ve P o si t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.13: MC Simulation, NB, Mean Relative Position Error Monte Carlo: Mean Relative Position Error with 1-sigma error bars: Noise & Bias, Random Initial Guess Error Case 1 10 100 1000 10000 100000 1000000 10000000 100000000 1357101520253035374045505560657075808590 Tracking Arc (days) M ean R e l a t i ve P o si t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.14: MC Simulation, NB, Mean Relative Position Error with sigmas 137 Figure 6.14 shows that there is an unperceivable standard deviation on the relative cases, whereas, on the non-Relative cases this deviation goes from +/- a few km to a few hundred meters. Figure 6.15 is an error bar plot of the mean relative velocity error with indicating very small standard deviations on the Relative cases after day 20. Monte Carlo: Mean Relative Velocity Error w/ 1-sigma error bars: Noise & Bias, Random Initial Guess Error Case 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 10 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) Me a n R e la tiv e V e lo c i ty E r r o r ( m /s ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.15: MC Simulation, NB, Mean Relative Velocity Error with sigmas The raw data for the mean relative position and velocity with the (+/-) 1-sigma standard deviations are summarized in Table 6.4 below. Of note, DSN has standard deviations of the same magnitude as the measurement itself, with the sigma-error on the order of km by day 90. Currently, there is an inter-spacecraft accuracy requirement of about 300 m [8], which may be tighten to 20 m or less in the future [3]. These findings suggest that it would be impossible to meet this requirement without Relative data cases (D/V/R and D/R). DSN/VLBI could only achieve a 300 meter relative position error accuracy after a ~83 day tracking arc. 138 Table 6.4: Summary of Mean Relative Position & Velocity with Standard Deviations Mean Relative Position Error (m) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 4278.9 (+/- 2224) 3.468 (+/-0.113) ------ 3.469 (+/- 0.111) 20 1475 (+/- 820.29) 3.474 (+/-0.063) ------ 3.485 (+/- 0.081) 30 1209.9 (+/- 580.48) 3.497 (+/- 0.051) 1.24 x 10 7 (+/- 7.06 x 10 6 ) 3.510 (+/- 0.049) 60 527.01 (+/- 286.1) 3.629 (+/- 0.054) 2.15 x 10 5 (+/- 1.18 x 10 5 ) 3.631 (+/- 0.054) 90 148.49 (+/-89.641) 3.698 (+/- 0.049) 17996 (+/- 12623) 3.698 (+/- 0.049) Mean Relative Velocity Error (m/s) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 0.06618 (+/- 0.0376) 0.000299 (+/- 1.53 x 10 -4 ) ------ 0.001447 (+/- 9.15 x 10 -4 ) 20 0.001611 (+/- 7.62 x 10 -4 ) 2.55 x 10 -6 (+/- 7.21 x 10 -7 ) ------ 2.54 x 10 -6 (+/- 7.18 x 10 -7 ) 30 8.54 x 10 -4 (+/- 3.92 x 10 -4 ) 2.46 x 10 -6 (+/- 2.45 x 10 -7 ) 2.47 (+/- 1.304) 2.46 x 10 -6 (+/- 2.46 x 10 -7 ) 60 1.63 x 10 -4 (+/- 1.13 x 10 -4 ) 2.33 x 10 -6 (+/- 4.97 x 10 -8 ) 5.00 x 10 -2 (+/- 0.026) 2.33 x 10 -6 (+/- 4.97x 10 -8 ) 90 4.79 x 10 -5 (+/- 2.39 x 10 -5 ) 1.60 x 10 -6 (+/- 2.32 x 10 -8 ) 4.19 x 10 -3 (+/- 0.002) 1.60 x 10 -6 (+/- 2.31x 10 -8 ) Without biases, (denoted as "noB") and initial guess error there is a dramatic improvement in the already outperforming DSN/VLBI/Relative and DSN/Relative cases; there is an improvement in DSN as well. Just random noise is applied on the simulation in this section. Figure 6.16 is the absolute position error and shows that DSN gets as good as DSN/VLBI by day 65. The Relative case errors continue to decrease rather than remain constant as with biases. 139 Monte Carlo Mean Absolute Position Error: Noise, no Bias, no Initial Guess Error 1 10 100 1000 10000 100000 1000000 10000000 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M ean A b so l u t e P o s i t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.16: MC Simulation, noB, Mean Absolute Position Error Figure 6.17 shows a lot more correlation between the error bars of DSN and DSN/VLBI around day 45. DSN/VLBI and DSN/Relative overlies from days 10-20 and then DSN/Relative and DSN/VLBI/Relative overlap for the rest of the time. This is a strong indicator that VLBI is no longer effective once these two curves overlap. It can also be inferred that VLBI no longer becomes effective once DSN decreases below DSN/VLBI. Figure 6.18 shows the mean relative position error with standard deviations. VLBI appears to lose effectiveness by day 65 when the DSN curve becomes better than the DSN/VLBI curve. Figures 6.19 and 6.20 show the absolute and relative velocity with their standard deviations. 140 Monte Carlo Mean Absolute Position Error with 1-sigma standard deviation error bars: Noise, no Bias, no Initial Guess Error 0.1 1 10 100 1000 10000 100000 1000000 10000000 1 3 5 7 101520253035374045505560657075808590 Tracking Arc (days) M e a n A b s o l u te P o s i ti o n E r r o r (m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.17: MC Simulation, noB, Mean Absolute Position Error with sigmas Monte Carlo Mean Relative Position Error with 1-sigma standard deviation error bars: Noise, no Bias, no Initial Guess Error 0.01 0.1 1 10 100 1000 10000 100000 1000000 10000000 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M ean R e l a t i ve P o si t i o n E r r o r ( m ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.18: MC Simulation, noB, Mean Relative Position Error with sigmas 141 Monte Carlo Mean Absolute Velocity Error with 1-sigma Standard deviation Error bars: Noise, no Bias, no Initial Guess Error 0.0000001 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 1 3 5 7 101520253035374045505560657075808590 Tracking Arc (days) M e a n A b s o lu te V e lo c i ty E r r o r (m /s ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.19: MC Simulation, noB, Mean Absolute Velocity Error with sigmas Monte Carlo Mean Relative Velocity Error with 1-sigma standard deviation error bars: Noise, no Bias, no Initial Guess Error 0.00000001 0.0000001 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) M ean R e l a t i ve V e l o ci t y E r r o r ( m / s ) DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative Figure 6.20: MC Simulation, noB, Mean Relative Velocity Error with sigmas Table 6.5 is a summary of the mean absolute position and velocity for the Monte Carlo simulation with noise and no biases. Note that the Relative-data corrects to much larger than 3 over the entire tracking arc for the Mean absolute position and the errors are 142 smaller than those seen in Table 6.3. Table 6.6 is a summary of the mean relative position and velocity for the Monte Carlo simulation presented in this section. Table 6.5: Summary of Mean Absolute Position & Velocity with Standard Deviations Mean Absolute Position Errors (m) (+/- 1-? ) Day s DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 5076.3 (+/- 1155.9) 2246 (+/- 1349.7 ) ------ 4414000 (+/- 3467000) 20 1717 (+/- 773.1) 412.63 (+/- 336.21) 8.14 x 10 5 (+/- 9.85 x 10 5 ) 732.7 (+/-608.37 ) 30 1356.5 (+/- 664.78) 154.17 (+/- 120.81) 44549 (+/- 33813) 193.47 (+/- 138.71 ) 60 590.8 (+/- 287.68) 14.471 (+/- 12.021 ) 844.9 (+/- 471.73) 14.322 (+/- 12.103 ) 90 238.61 (+/- 110.54) 3.5992 (+/- 2.1706 ) 75.061 (+/- 57.118) 3.58 (+/- 2.1834 ) Mean Absolute Velocity Errors (m/s) (+/- 1-? ) Day s DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 0.065252 (+/- 0.024914 ) 0.003078 (+/- 0.002031) ------ 0.17008 (+/- 0.12875) 20 0.001674 (+/- 0.000748 ) 3.24 x 10 -5 (+/- 2.52 x 10 -5 ) 0.1521 (+/-0.15522) 5.66 x 10 -5 (+/- 4.59 x 10 -5 ) 30 0.000812 (+/- 0.000364 ) 1.58 x 10 -5 (+/- 1.06 x 10 -5 ) 0.009395 (+/- 0.006804) 1.84 x 10 -5 (+/-1.33 x 10 -5 ) 60 0.000175 (+/- 8.01 x 10 -5 ) 2.48 x 10 -6 (+/- 1.69 x 10 -6 ) 0.000195 (+/- 0.000117) 2.48 x 10 -6 (+/- 1.69 x 10 -6 ) 90 4.49 x 10 -5 (+/- 2.26 x 10 -5 ) 8.25E-07 (+/- 5.20E-07) 1.47 x 10 -5 (+/- 1.00 x 10 -5 ) 8.23 x 10 -7 (+/- 5.23 x 10 -7 ) 143 Table 6.6: Summary of Mean Relative Position & Velocity with Standard Deviations Mean Relative Position Errors (m) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 3595.6 (+/- 1430.6 ) 0.16899 (+/- 0.0584 ) ------ 0.17643 (+/- 0.055822) 20 1449.1 (+/- 767.26) 0.096109 (+/-0.043723 ) 6.31 x 10 5 (+/-8.54 x 10 5 ) 0.10383 (+/- 0.050238) 30 1197.8 (+/- 659.26 ) 0.078791 (+/0.027683 ) 36700 (+/-27589 ) 0.081075 (+/- 0.02763 ) 60 503.98 (+/- 327.69 ) 0.065993 (+/- 0.028862) 699.78 (+/- 487.49) 0.065739 (+/- 0.029286) 90 239.51 (+/- 150.56) 0.057915 (+/- 0.026587) 76.138 (+/- 59.561 ) 0.057976 (+/- 0.02658) Mean Relative Velocity Errors (m/s) (+/- 1-? ) days DSN/VLBI DSN/VLBI/Relative DSN DSN/Relative 1 0.062469 (+/- 0.032609 ) 0.00029989 (+/- 0.00014873 ) ------ 0.001255 (+/- 0.000961) 20 0.001732 (+/- 0.00102 ) 9.51x 10 -7 (+/-5.21 x 10 -7 ) 0.14373 (+/- 0.16548) 9.47x 10 -7 (+/-5.18 x 10 -7 ) 30 0.000816 (+/- 0.000461) 3.49 x 10 -7 (+/-1.67 x 10 -7 ) 0.009325 (+/- 0.007134) 3.49 x 10 -7 (+/- 1.67 x 10 -7 ) 60 0.00018 (+/- 9.97 x 10 -5 ) 7.80 x 10 -8 (+/- 3.69 x 10 -8 ) 0.000195 (+/- 0.000121) 7.81 x 10 -8 (+/- 3.69 x 10 -8 ) 90 4.81 x 10 -5 (+/- 2.82 x 10 -5 ) 3.32 x 10 -8 (+/- 1.47 x 10 -8 ) 1.59 x 10 -5 (+/- 1.33 x 10 -5 ) 3.32 x 10 -8 (+/- 1.47 x 10 -8 ) The next two figures predict the effect that increasing the number of measurements per day from 1 to 3 would have on each case. A Monte Carlo dataset was collected for 1 and 3 measurements per day at day 20 with noise and bias and no initial guess errors to determine if more frequent measurements better the OD. Figure 6.21 shows the absolute position error, indicating an error reduction by hundreds of meters for most cases. The least improvement is seen on the DSN/VLBI/Relative case. The greatest improvement is on DSN. The standard deviation bars are also plotted on the columns and the standard 144 deviation remains comparable for all the cases except DSN/VLBI/Relative in which it decreases for three measurements per day. Frequency of Measurements Comparison: Absolute Position Error Monte Carlo noise, Bias, no initial guess error at 20 days 1820.17 1141.44 383.30 296.24 1548668.06 832009.50 900.04 658.71 1 10 100 1000 10000 100000 1000000 10000000 M ean A b so l u t e P o si t i o n E r r o r w i t h S t an dar d D e vi at i o n s (m ) l o g - y D/V_1 D/V_3 D/V/R_1 D/V/R_3 D_1 D_3 D/R_1 D/R_3 Figure 6.21: Frequency of Measurements, Absolute Position Error In Figure 6.22, increasing the measurements per day for mean relative position error has very little effect on the Relative cases, as both bars are nearly equal in value for DSN/VLBI/Relative and DSN/Relative. There is also insignificant standard deviation in both columns. However, for the non-Relative cases, there is a decrease in error with the increased measurement availability. 145 Frequency of Measurements Comparison: Relative Position Error Monte Carlo noise, Bias, no initial guess error at 20 days 1820.17 1043.24 3.465 3.462 1207920.64 643893.77 3.474 3.478 1 10 100 1000 10000 100000 1000000 10000000 M ean R e l a t i ve P o si t i on E r r o r w i t h S t an dar d D e vi at i o n s (m ) l o g - y D/V_1 D/V_3 D/V/R_1 D/V/R_3 D_1 D_3 D/R_1 D/R_3 Figure 6.22: Frequency of Measurements, Relative Position Error 6.4. Discussion of Results Summary comparisons of the Monte Carlo runs are exhibited in the following bar graphs. Additional comments will be made concerning various requirement and constraints. Figure 6.23 illustrates how many days it takes each case to get to an absolute position accuracy of 1 kilometer for both Monte Carlo runs. Linear interpolation was used to determine exactly when each case would reach 1 kilometer. Without biases ("noB"), all cases reduce to 1 kilometer error sooner than with biases ("NB"). There is only a slight improvement for DSN/VLBI and the Relative cases. DSN however, improves greatly, though this can be attributed to no initial guess error since DSN was shown to handle biases quite well. Note that for the Monte Carlo Noise/Bias Run, the DSN-only case never achieved an error as low as 1 kilometer. In order to obtain the projected 99.63 day 146 bar, linear interpolation was applied to best estimate when this would be achieved beyond the 90 day commissioning period. Monte Carlo: Days Until Reach 1 Kilometer in Absolute Position Error 36.11 36.84 6.82 7.35 57.76 99.63 19.10 19.17 0 20 40 60 80 100 120 D ays DSN/VLBI: noB DSN/VLBI: NB DSN/VLBI/Relative: noB DSN/VLBI/Relative: NB DSN: noB DSN: NB DSN/Relative: noB DSN/Relative: NB Figure 6.23: Summary Chart of Days Until Absolute Position Error Equals 1 Kilometer For relative position error, only the Relative cases (DSN/VLBI/Relative and DSN/Relative) are capable of achieving an accuracy of < 20 meters in the 90 day time period. Figure 6.24 illustrates how many days it takes to obtain a relative position accuracy of the less stringent 300 meters for both Monte Carlo runs. Again linear interpolation was used and some values are projected beyond the 90 day period. Within 1 day, both DSN/VLBI/Relative and DSN/Relative are less than 300 meters. 147 Monte Carlo: Days Until Reach 300 Meters in Relative Position Error 77.27 81.69 1.00 1.00 69.81 101.18 1.00 1.00 0 20 40 60 80 100 120 D ays DSN/VLBI: noB DSN/VLBI: NB DSN/VLBI/Relative: noB DSN/VLBI/Relative: NB DSN: noB DSN: NB DSN/Relative: noB DSN/Relative: NB Figure 6.24: Summary Chart of Days Until Relative Position Error Equals 300 meters A current requirement is that the LISA spacecraft radius, or semi-major axis should be between 0.9 AU ( = 1.346 x 10 11 m) and 1.1 AU (= 1.646 x 10 11 m) [8]. This constraint is met as each LISA for all cases stays within 2.83 x 10 8 m of the nominal semi-major axis. There is currently a more stringent requirement on the period, such that each LISA should be within 38 second of nominal heliocentric period before beginning drag free science operations [8]. Since the period is related to the semi-major axis (i.e. P = 2?* ? 3 a ), the semi-major axis must be controlled to < 60 kilometers of the desired. The prior stated 2.83 x 10 8 m difference from the nominal semi-major axis translates to a 1.04 day maximum period difference. Since LISA?s orbital period is one year, this period difference is merely 0.284% of the nominal period. Curiously, the earlier stated constraint of the semi-major axis being between 0.9 AU and 1.1 AU translates into a much larger period offset of 53-56 days of the nominal, or a 14.6%-15.4% error. This 148 simulation achieved semi-major axis within the loose constraint, while the periods were different by up to one day. Future work must investigate closing the gap on the tight period requirement through the use of ?V maneuvers prior to dropping the propulsion modules. Parameter uncertainties decrease with increasing time as is evidenced by Figure 6.25, which shows the parameter uncertainty of the DSN/VLBI/Relative case for all three LISA spacecraft and all three position components. The Z-component, or the cross-track axis of the satellite RSW system for all three spacecraft, remains the most uncertain (the top three lines of the Figure 6.25) in low declination orbits. The Y- or along-track components are the middle three lines, while the X- or radial components are the bottom three lines. A longer tracking arc ensures more certainty in all components. DSN/VLBI/Relative Mean Parameter Uncertainty for All LISA Position Components 0 100 200 300 400 500 600 1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90 Tracking Arc (days) Pos i t i on C o m pone n t s ( m ) X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 Figure 6.25: MC (NB) Parameter Uncertainty of LISA 1 position-components 149 The parameter uncertainties (PU) and the mean estimated state ( ) ? (Xmean ) at day = 7 for the Monte Carlo (Noise, Bias) case are listed in Table 6.7. Table 6.7: Parameter Uncertainties (PU) for DSN/VLBI/Relative Case at day 7 LISA 1 LISA 2 LISA 3 ) ? ( 1 Xmean +/- PU 1 ) ? ( 2 Xmean +/- PU 2 ) ? ( 3 Xmean +/- PU 3 X (m) 5.440 32.411 2.488 x 10 9 47.786 -2.488 x 10 9 35.773 Y (m) 1.481 x 10 11 198.210 1.503 x 10 11 260.680 1.503 x 10 11 177.120 Z (m) 2.476 x 10 9 286.950 -1.292 x 10 9 381.750 -1.292 x 10 9 261.200 X & (m/s) -3.007 x 10 4 3.592 x 10 -5 -2.964 x 10 4 4.701 x 10 -5 -2.964 x 10 4 3.327 x 10 -5 Y & (m/s) 3.333 x 10 -6 2.194 x 10 -5 2.476 x 10 2 3.170 x 10 -5 -2.476 x 10 2 2.269 x 10 -5 Z & (m/s) 3.624 x 10 -6 3.965 x 10 -5 4.269 x 10 2 5.454 x 10 -5 -4.269 x 10 2 3.851 x 10 -5 150 Chapter 7: Conclusion This thesis has studied the Orbit Determination (OD) problem for the proposed Laser Interferometer Space Antenna (LISA) mission in which three satellites in heliocentric orbit must maintain a triangular formation. This work investigated the impact of augmenting traditional Earth-based DSN and VLBI OD data with inter-spacecraft range and range rate (also termed "Relative" data) measurements available for the LISA mission. A weighted batch least squares algorithm was used to estimate spacecraft orbital parameters from four test combinations of measurement types: 1) DSN and VLBI, 2) DSN, VLBI and Relative, 3) DSN only, and 4) DSN and Relative. An extensive set of simulation-based experiments were conducted to study OD accuracy statistics. The two main modes of analysis to compare the test cases were the Covariance and Monte Carlo methods. Very Long Baseline Interferometry (VLBI) has been shown to improve initial data quality, particularly when DSN is the only other measurement available. However, when applied with a noise of 50 nrad, VLBI becomes obsolete after a certain amount of time since long arcs of range and Doppler data can provide angular position accuracy of 40 nrad or better [41]. Though VLBI or some type of angle data is helpful initially, LISA passes through zero declination and has periods of poor observability over the course of its orbit. LISA needs a longer Doppler arc, or tracking must be supplemented with an 151 alternative technique for measuring angles to become observable. VLBI gives orientation to the range measurement when Doppler is not good initially. More measurements per day were found to improve the absolute position accuracy. Given LISA?s distance from the Earth, its orbital period of one year, and LISA unique geometry with respect to Earth, a sufficient tracking arc length tells more about the spacecraft than denser data. Future work would be needed to investigate the concept of more frequent data further. Confidence in the answer also improves with longer tracking arcs. Parameter uncertainties and standard deviations in the covariance matrix diagonals associated with each element of position and velocity as well as standard deviations associated with measured errors (i.e. absolute position and relative error) all decrease with an increasing tracking arc. The most major finding of this work is that Relative data is worth including to improve not only the relative knowledge of the spacecraft, but also for absolute position and velocity OD as well. There was a greater than expected 3 improvement seen for absolute OD accuracy for the cases with Relative data. VLBI is also worth incorporating as it does a better job initially of tracking LISA in its low declination orbit. After about 20-35 days, VLBI can be reallocated to other assets, and the tracking can rely on DSN and Relative data without loss of accuracy. DSN-only needs a 20 day tracking arc with a good initial guess; if the initial guess is off by more than 1 km, a 30-35 tracking arc is needed. The other cases could correct for initial guess errors up to 10 km, though no check was done on the upper limit. Noises in the measurements are filtered out with a 152 sufficiently long tracking arc, but the biases are constant. The biases affected the Relative data more than the non-Relative, which points to a need for the calibration and the subtraction of known biases for better results. The current inter-spacecraft accuracy requirement of 300 m may be tightened in the future to an accuracy of 20 m or less. DSN/VLBI can barely make that requirement until 85 days of tracking. If the requirement were more stringent, then only the Relative are capable of achieving that level of accuracy. The semi-major axis constraint is met by all cases, but the more stringent 38 second period error constraint prior to propulsion module drop-off is still harder to achieve. This simulation estimates error to be within one day of the nominal period. The overall best case from this study is DSN/VLBI/Relative, i.e. using all available measurements, including the very crucial inter-spacecraft range and range rate. After 7 days of tracking LISA, DSN/VLBI/Relative is already achieves below one kilometer absolute OD accuracy and less than 3.5 meters in relative OD error. 7.1. Future Work Future work would include simulating a more realistic model of LISA. The perturbations of planets would need to be added for a more complete and accurate force model of LISA in its space environment. The nominal orbit of Dhurandhar [18], who did substantial 153 work with the Clohessy-Wiltshire equations on the unique +/- 60 degree tilt of the LISA orbit in a perturbed environment, would most likely be used. The rotation of Earth and the tracking station location on Earth's surface would also have to be modeled. In particular, case studies are required that use the location of specific tracking stations, their visibility of the spacecraft, and sensors that will actually be used. Even more specific simulations would include trade studies to determine which sensors are best, the choice of a Deep Space Station antenna(s) to be used, VLBI antennas, and the selection of the quasar to be used as a reference for angular measurements. More investigation over the ideal amount of measurements per day could be conducted. It as found that for this orbit increasing the number of measurements to three per day does improve the absolute OD accuracy and to a lesser degree the relative OD at 20 days. However, could the accuracy stay the same for skipping days? For example, consider a three-day rotation where LISA 1 is tracked once per day on day 1, LISA 2 on day 2, etc. The data combination of DSN/VLBI/Relative was found to be the best; however, angle data becomes ineffective after 20-35 days. Further combinations that switch data types on and off when they are or are not effective could be investigated to see if the results are as accurate as the DSN/VLBI/Relative case. For example, using DSN/VLBI/Relative at the beginning and then turning off VLBI to just rely on DSN/Relative when the angles are no longer effective would be an intuitive choice. Cost related studies are also a spin-off of this thesis. The practical implication of using different measurements as a function of cost can be studied, a total of six measurement 154 types (range, range rate, right ascension, declination, inter-spacecraft range and range rate) were used here. The potential cost savings of turning VLBI off when no longer effective could be investigated, since this data type requires two antennas on Earth to obtain its measurements. The possibility of using less DSN and relying more on VLBI and inter-spacecraft data is also of interest, since the current configuration of DSN is struggling to support the deep space communications demands of the coming decades. However, VLBI could potentially pose an advantage over DSN since it uses two smaller antennas rather than a larger single antenna, which saves in operating and production costs [42]. Further work with weighting the data, using less or more noise and different biases also could help understand other nuances of this problem, the measurements, and how the error sources interact. Biases were found to have a greater effect on the relative data, therefore, estimating biases and working to better calibrate these biases in the instruments could improve OD. The 90 day time period of the close to nominal commissioning period was used as the thesis benchmark. Future work could also include looking at other specific OD events around the three month time period such as when the propulsion module drop-off and will "kick" the spacecraft, and how long it will take to return back to the nominal orbit after that disturbance. The ideal tracking arcs should be determined before and after propulsion module drop-off. Lastly, the use of Relative i.e. inter-spacecraft data types could be very beneficial to future constellation missions as a method to get increased OD 155 accuracy, which is an improvement that cannot be done by Earth-based tracking alone. Better OD accuracy and navigation means better science and mission success. 156 Appendices Appendix A: Supporting Calculations Newton's Law of Gravitation: () 1 22 21 21 3 21 21 111 , GMwhere rr mGm RR RR mGm RmF = ? = ? =? ? ? == ? ? && and the mass, m 2 of any LISA spacecraft is negligible compared to the larger mass (i.e. Earth, Jupiter and the Sun in this case). The laws of sines and cosines were used to find the distances between the spacecraft and the planetary body. Law of sines: () () ()??? sinsinsin cba == Law of cosines: a 2 = b 2 + c 2 ? 2*b*c*cos(?) ? Force of Sun on LISA: 2 6 2 2311 2 1093.5 )870,597,149( /1032712428.1 s km km skm r F SL S ? ??? ?? = ? = ? where, the distance of LISA to the Sun is: r SL = 1 Astronomical Unit = 1 AU = 149,597,870 km ? Force of Earth on LISA: 2 10 2 235 2 106555.1 )6708.043,068,49( /10986004415.3 s km km skm r F EL E ? ??? ?? = ? = ? where, the distance of the closest LISA from the Earth is: r EL = 49,068,043.6708 km ? Force of Jupiter on LISA: 2 10 26 238 2 106735.3 )346.15248462110740( /10268.1 s km kmkm skm r F JL J ? ??? ?? ?? = ? = ? where Jupiter at closest approach to Sun is 740 million km, and the furthest configuration of LISA from the Sun is 152,484,621.346 km. Compared to the Sun, the force of Earth on LISA is 0.00279 % and the force of Jupiter on LISA is 0.006%. 157 Appendix B: MATLAB CODES Table of Contents Orb2X2.m ?????????????????????????????158 -This script transforms a set of Keplerian orbital elements to Cartesian format. X2Orb.m ??..???????????????????????????160 -This script transforms a set of Cartesian elements to Keplerian orbital elements. LISA_EOMderivs.m ????????????????????????..164 -Propagates the input state using the equations of motion. NomOrb.m ????????????????????????????. 164 -Determines a nominal trajectory from Steve Hughes initial states. This script generates the ephemeris for one year by using the equations of motion. LISA_STMderivs.m ????????????????????????...166 -Maps the state transition matrix to each time in the propagation H_partials.m ?????????????????..?????????? 167 -This program presents the partials of the observation-state mapping matrix for each element for each state component. RealObs.m ????????????????????????????...171 -Calculates the Earth-based observations (range, range rate, right ascension and declination) at each measurement time, it also simulates the "observed" measurements by including noises and biases. InterSCD.m ????????????????????????????174 -Calculates the inter-spacecraft observations (inter-spacecraft range and range rate) at each measurement time, it also simulates the "observed" measurements by including noises and biases. LISA_WBLS.m ??????????????????????????. 176 -This program is the weighted batch least squares algorithm for each of the four cases. MCRUNbp.m ??????????????????????????? 188 -This program calls LISA_WBLS.m and runs the Monte Carlo simulation and computes relevant statistics. 158 Appendix B: MATLAB CODES Orb2X2.m (from Section 5.2.2.The Nominal Orbit) function [R, V] = Orb2X2(Ephem, t, mu, type) %*********************************************************** % %INPUT: Ephem1 = [a e i Om w tau], the six classical orbital elements % [semi-major axis, eccentricity, inclination, longitude of the % ascending node, argument of periapsis, time of periapsis] % % Ephem2 = [a e i Om w nu], the six classical orbital elements % [semi-major axis, eccentricity, inclination, longitude of the % ascending node, argument of periapsis, true anomaly] % % Ephem3 = [a e i Om ~omega L] % ~omega = longitude of perhelion % L = mean longitude % % Ephem4 = [a e i Om w M], the six classical orbital elements % [semi-major axis, eccentricity, inclination, longitude of the % ascending node, argument of periapsis, mean anomaly] % % t = time of epoch associated with the particular ephemeris % mu = gravitational parameter % type = 1 = Program uses Ephem1 as input % = 2 = Program uses Ephem2 % = 3 = Program uses Ephem3 % = 4 = Program uses Ephem4 % %OUTPUT : R - Position Vector % V - Velocity Vector % %CALLS : newton_meth.m, succ_sub_hyperbola.m % %Ephem2 = [149578965033 .01034342962 1.023240970 212.7394780 %90.76228189 136.5015274];[R, V] = Orb2X2(Ephem2, %2.456262500000000e+006, 389600.4415, 2) % %mu =4*(pi^2);t=.010576712; [R, V] = Orb2X(Ephem,t, mu) %REFERENCE: Equations from Battin(pp. 123-125), Tolson, Chobotov %Ephem = [6828973.232519 .0090173388450585 28.474011884869 %35.911822759495 -44.55584705279 43.8860381032208] %*********************************************************** %***INITIALIZATION***************************************** if nargin < 4 type = 1 %assume type 1 ephemeris inputs end if nargin < 3 mu = 3.9860044e5 %km^3/s^2, assume gravitational parameter of the Earth end r2d = 180/pi; d2r = pi/180; %*********************************************************** a = Ephem(1); %semi-major axis (usu. m/s or km/s) if a > 0 disp('elliptic orbit') elseif a < 0 disp('hyperbolic orbit') else %z = 1/a = 0 disp('parabolic orbit') end e = Ephem(2); %eccentricity (unitless) i = Ephem(3)*d2r; %inclination, from deg. to radians Om = Ephem(4)*d2r; %Right Ascension of the Ascending Node, from deg. %to radians if type ~= 3 %thus type 1, 2, or 4 w = Ephem(5)*d2r; %argument of periapsis, from deg. to radians w_bar = Om + w; %longitude of periapsis = argument of periapsis - %longitude of the ascending node (radians) 159 elseif type == 3 w_bar = Ephem(5)*d2r; %longitude of perihelion w = w_bar - Om; %argument of periapsis end if e == 0 %circular orbit, no periapsis w_bar = NaN; %undefined end if e < 1 mean_mo = sqrt(mu/a^3); %mean motion elseif e > 1 mean_mo = sqrt((-1*mu)/a^3); %mean motion else %e == 1 p = a*(1-e^2); %semi-latus rectum mean_mo = sqrt(mu/p^3); %mean motion end p = a*(1-e^2); %semi-latus rectum %***Solve Kepler's equation for the Eccentric anomalies for type 1 and 3 tol=1e-10; if type == 1 tau = Ephem(6); %This Ephemeris uses the time of perapsis passage as %the 6th element %Find Mean Anomaly ep = w_bar - (mean_mo*tau); %the mean longitude at epoch time t=0 M = (mean_mo * t) + ep - w_bar; else type == 3 L = Ephem(6)*d2r; M = L - w_bar; tau = t - (M/mean_mo); if e < 1 %ellipse, solve Kepler's equation for E, Eccentric Anomaly %Newton's iteration method. Call to file "newton_meth.m" [E,its] = newton_meth(e,M,tol); E nu = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); %radians elseif e > 1 %hyperbola, solve Kepler's equation for a hyperbola for F %Successive Subsitution method, call "Succ_sub_hyperbola.m" [F,its] = succ_sub_hyperbola; F else %e == 1 %parabola P = 0 %NO Eccentric Anomaly for Parabolic case end end if type == 4 M = Ephem(6)*d2r tau = t - (M/mean_mo) if e < 1 %ellipse, solve Kepler's equation for E, Eccentric Anomaly %Newton's iteration method. Call to file "newton_meth.m" [E,its] = newton_meth(e,M,tol); E nu = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); %radians elseif e > 1 %hyperbola, solve Kepler's equation for a hyperbola for F %Successive Subsitution method, call "Succ_sub_hyperbola.m" [F,its] = succ_sub_hyperbola; F else %e == 1 %parabola P = 0 %NO Eccentric Anomaly for Parabolic case end end %Solve for true anomaly %if type == 1|3 % nu = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); %radians %end %***Find the ECCENTRIC and MEAN ANOMALIES for type 2 160 if type == 2 nu = Ephem(6)*d2r; %Otherwise the 6th element is the true anomaly, convert input into radians if e < 1 cE = (e + cos(nu))/(1+(e*cos(nu))); sE = (sqrt(1-e^2)*sin(nu))/(1+e*cos(nu)); E = atan2(sE,cE) M = E - e*sin(E); elseif e > 1 cF = (e + cos(nu))/(1+(e*cos(nu))); sF = (sqrt(e^2-1)*sin(nu))/(1+e*cos(nu)); F = atan2(sF,cF) M = e*sinh(F)-F; else %e==1 P = 0; M = (1/2)*tan(nu/2) + (1/6)*(tan(nu/2)^3); end end %ML = w_bar + M; %mean longitude %M = ML - L; %Mean anomaly %***CALCULATE POSITION (R) & VELOCITY (V)************ if e < 1 rmag = a*(1-e*cos(E)); %magnitude position for an ellipse elseif e > 1 rmag = a*(1-e*cosh(F)); %magnitude position for a hyperbola else % e==0 rmag = (p/2)*(sec(nu/2)^2);%magnitude position for a parabola end theta = (w + nu); %argument of latitude = (w + nu) in radians sN = sin(Om); cN = cos(Om); sT = sin(theta); cT = cos(theta); si = sin(i); ci = cos(i); sw = sin(w); cw = cos(w); %position vector R = [rmag*((cT*cN)-(sT*sN*ci)); rmag*((cT*sN)+(sT*cN*ci)); rmag*(sT*si)]; h = sqrt(mu*p); %angular momentum scalar mu_h = mu/h; %velocity vector V = [-mu_h*((cN*(sT + e*sw)) + (sN*(cT + e*cw)*ci));... (-mu_h*((sN*(sT + e*sw))-(cN*(cT + e*cw)*ci))); (mu_h*((cT + e*cw)*si))]; X2Orb.m (from Section 5.2.2. The Nominal Orbit) function [E] = X2Orb(R, V, t, mu) %*********************************************************** %function [Ephem, Ephem_2, Misc, E] = X2Orb(R, V, t, mu) % %Function to obtain the six orbital elements plus some other orbital %information from an inputted position and velocity vector. % %INPUTS : R = Position Vector (m) % V = Velocity Vector (m/sec) % t = epoch time of inputted R and V vector (sec) % mu = inputted gravitational parameter (m/sec^3) % %OUTPUTS : Can customize to read out different sets of elements. % Ephem = [a e i Om w tau], the six classical orbital elements % [semi-major axis, eccentricity, inclination, longitude of the % ascending node, argument of periapsis, time of periapsis] in 161 % deg, time and length units % Ephem_2 = [f(True anomaly), M(Mean anomaly), ... % w_bar_true(true longitude of periapsis),... % w_bar(longitude of periapsis), ... % uo(argument of latitude at epoch), lo(true longitude at epoch)] % Misc = [h(angular momentum), p(semi-latus rectum), ra(apoapsis), %rp(periapsis),... % Period], extra quantities to display % E = [a e i w Om M f Period] % %REFERENCE: Equations from (1)"Fundamental of Astrodynamics" by %Bate, Mueller % & White and Battin; (2)"Statistical Orbit Determination" by Byron Tapley et. % al.; (3)"Orbital Mechanics" Chobotov (4) Vallado, ELORB algorithm % pp.146-147 %*********************************************************** %INITIALIZATION if nargin < 4 mu = (1.32712428e20); %m3/s2, Sun's gravitational parameter; %1.327124e11*(1000^3); %m3/s2, Sun's gravitational parameter %mu = 3.9860044e5*(1000^3) %m^3/s^2, assume gravitational %parameter of the Earth %error('Must pass in at least R,V, t and mu. Type "help X2Orb" for more %information.'); end rmag = norm(R); %magnitude of position vector vmag = norm(V); %magnitude of velocity vector r2d = 180/pi; %conversion from radians to degrees d2r = pi/180; %conversion from degrees to radians %ANGULAR MOMENTUM VECTOR, hvec hvec = cross(R,V); %hvec is perpendicular to the plane of the orbit hx = (R(2)*V(3)) - (V(2)*R(3)); %the x,y,z,xy components of the h vector hy = (R(3)*V(1)) - (V(3)*R(1)); hz = (R(1)*V(2)) - (V(1)*R(2)); hxy = (hx^2 + hy^2)^(1/2); h = norm(hvec); %the scalar of the h vector %NODE VECTOR, nvec, pointing along the line of nodes in the direction %of the ascending node kvec =[0 0 1]; %k-vector nvec = cross(kvec, hvec); %node vector n = norm(nvec); %node W = cross(R,V)./norm(cross(R,V)); W_check = hvec./norm(hvec); %ECCENTRICITY VECTOR, evec, points from focus of orbit to the %periapsis evec = (1/mu)*(((vmag^2 - (mu/rmag))*R) - ((dot(R,V))*V)); e = norm(evec); Vs = dot(V,V); vs2 = vmag^2; Energy = (Vs/2) - (mu/rmag) ; %Vis-Viva Equation evec_check = cross(V,(hvec./mu)) - (R./rmag); e_check = (1 + ((2*Energy*h^2)/mu^2))^(1/2); %SEMI-LATUS RECTUM & SEMI-MAJOR AXIS if e~= 0 a = -mu/(2*Energy); %semi-major axis p = a*(1-e^2); %semi-latus rectum rp = p/(1+e); %periapsis ra = p/(1-e); %apoapsis a_check = (rp + ra)/2; a_check2 = 1/((2/rmag)-(vmag^2/mu)); else 162 a = inf; % semi-major axis equals infinity p = h^2/mu;%semi-latus rectum end %PERIOD Period = ((2*pi)*(a^(3/2)))/(sqrt(mu)); %seconds Period_check = ((2*pi)*(a_check^(3/2)))/(sqrt(mu)); %INCLINATION, i i = acos(dot(hvec,kvec)/h)*r2d; %all 3 methods should yield same value i_check = acos(hz/h)*r2d; i_check2 = acos(dot(W,kvec))*r2d; if i <= 90; disp('direct, easterly orbit, inclination between 0 and 90 deg'); elseif i <= 180; disp('retrograde orbit, inclination is between 90 and 180 deg.'); else %i>180, %inclination is always less than 180 deg. disp('error, i is always less than 180 deg.'); end %LONGITUDE OF THE ASCENDING NODE, Om ivec = [1 0 0]; N = (cross(kvec,W))./norm(cross(kvec,W)); nvec(2); if nvec(2) < 0 %if the jth element of the node vector is < 0 then Om = 360deg - Om Om_pre = acos(nvec(1)/n)*r2d; Om = 360 - Om_pre; else Om = acos(nvec(1)/n)*r2d; end %TIME OF PERIAPSIS z = 1/a; if z > 0 %elliptical case mean_mo = sqrt(mu*z^3); cE = inv(e)*(1-z*rmag); sE = z^2*(dot(R,V))/(mean_mo*e); E = atan2(sE,cE); M = (E - (e*sin(E)))*r2d; M_check = (E - (e_check*sin(E)))*r2d; tau = t -(E - e*sin(E))/mean_mo; tau_check = t - (M_check*d2r)/mean_mo; elseif z < 0, %hyperbolic case mean_mo = sqrt(mu*-z^3); cE = inv(e)*(1-z*rmag); sE = z^2*(dot(R,V))/(mean_mo*e); E1 = asinh(sE); E2 = acosh(cE); E = E2; if E1 < 0, E = -E2; end tau = t -(e*sinh(E) - E)/mean_mo; else % z == 0 %parabolic case mean_mo = sqrt(mu/p^3); cE = p/rmag - 1; sE = (dot(R,V))/rmag*sqrt(p/mu); E = atan2(sE,cE); tau = t -((tan(E/2)^3)/6 +tan(E/2)/2)/mean_mo; end %M %The following 3 elements (nu, uo, lo) can be substituted for tau %TRUE ANOMALY, the angle in satellite's orbit plane between periapsis %and satellite's position at a particular epoch, to if dot(R,V) < 0, %if the dot product of R and V is greater than 0, then the true anomaly should be less than 180 deg. fpre = acos((dot(evec,R))/(e*rmag)) *r2d; f = 360 - fpre; else 163 f = acos((dot(evec,R))/(e*rmag)) *r2d; end %ARGUMENT OF PERIAPSIS, w if e == 0, %circular orbit, no periapsis w = undefined; end if evec(3) < 0, %if the kth element of the eccentricity vector is >0 then w is %<180 deg. wpre = acos((dot(nvec,evec))/(n*e))*r2d; w = 360 - wpre; else w = acos((dot(nvec,evec))/(n*e))*r2d; end %LONGITUDE OF PERIAPSIS, w_bar, can be substituted for w w_bar = Om + w; if e == 0, %circular orbit, no periapsis w_bar = undefined; end %TRUE LONGITUDE OF PERIAPSIS, w_bar_true if evec(2) < 0 w_bar_true_pre = acos(evec(1)/e)*r2d; w_bar_true = 360 - w_bar_true_pre; else w_bar_true = acos(evec(1)/e)*r2d; end %ARGUMENT OF LATITUDE, uo, the angle between nvec and R at %epoch, to if R(3) < 0, %if the kth element of the position vector is < 0, then uo is < 180 deg. uo_pre = acos((dot(nvec,R))/(n*rmag)) *r2d; % = w + nu uo = 360 - uo_pre; else uo = acos((dot(nvec,R))/(n*rmag)) *r2d ;% = w + nu end %TRUE LONGITUDE, lo, the angle between I and ro(the radius vector to %the satellite at epoch, to), measured eastward to the ascending node if it %exists and then in the orbital plane to ro. If Om, w, and f are all %defined then lo. if Om == NaN, %equatorial orbit lo = w_bar + f; elseif w == NaN, %circular orbit lo = Om + uo; else lo = acos(R(1)/rmag)*r2d; if R(2) < 0 lo_pre = acos(R(1)/rmag)*r2d; lo = 360 - lo_pre; end end format long e %***OPTIONAL EPHMERIDES TO DISPLAY********************* %disp('Ephem = [a e i Om w tau]'); %Ephem = [a e i Om w tau]; %disp('Ephem_2 = [f(True anomaly), M(Mean anomaly), w_bar_true(true %longitude of periapsis), w_bar(longitude of periapsis),uo(argument of %latitude at epoch), lo(true longitude at epoch)]'); %Ephem_2 = [f M w_bar_true w_bar uo lo]; %disp('Misc =[h(angular momentum), p(semi-latus rectum), ra(apoapsis), %rp(periapsis),Period]') %Misc = [h, p, ra, rp, Period]; disp('Ephemeris = [a e i Om w M f Period]') E = [a e i Om w M f Period]; 164 LISA_EOMderivs.m (from Section 5.2.2. The Nominal Orbit) function Xdot = LISA_EOMderivs(t, X) %*********************************************************** %Xdot = LISA_EOMderivs(t, X) % %Simulate the orbits of 3 sun-orbiting LISA spacecraft with the given intial conditions %using the Equation of Motion(EOM). EOM = -(mue/(norm(r))^3).*r % %INPUT: X = initial state of the orbit = % [S/C #1 position & velocity, S/C #2 position & velocity, S/C #3 position & velocity, % Tracking Station Position & Velocity] in 3-D (1 x 24) % t = initial time % %OUTPUT: Xdot = the derivative of the inputted state = % Xdot = [rdot_1; rddot_1; rdot_2; rddot_2; rdot_3; rddot_3; rdot_rts; %rddot_rts] % %Command Window: [t, X] =ode45('LISA_EOMderivs', timespan, XO, %odeset('RelTol', 1e-12,'AbsTol',1e-12)); %*********************************************************** mus = (1.32712428e20); %m3/s2, Sun's gravitational parameter r1 = X(1:3); %m, LISA spacecraft #1 location rmag_1 = norm(r1); rdot_1 = X(4:6); %m/sec, LISA spacecraft #1 velocity rddot_1 = -(mus/(rmag_1^3))*r1; %m/s2, EOM for spacecraft #1 r2 = X(7:9); %m, LISA spacecraft #2 location rmag_2 = norm(r2); rdot_2 = X(10:12); %m/sec, LISA spacecraft #2 velocity rddot_2 = -(mus/(rmag_2^3))*r2; %m/s2, EOM for spacecraft #2 r3 = X(13:15); %m, LISA spacecraft #3 location rmag_3 = norm(r3); rdot_3 = X(16:18); %m/sec, LISA spacecraft #3 velocity rddot_3 = -(mus/(rmag_3^3))*r3; %m/s2, EOM for spacecraft #3 rts = X(19:21); %m, tracking station location IS Earth (treat Earth as a %point) rts_mag = norm(rts); rdot_rts = X(22:24); %m/sec, Earth velocity rddot_rts = -(mus/(rts_mag^3))*rts; %m/s2, EOM for spacecraft #3 Xdot = [rdot_1; rddot_1; rdot_2; rddot_2; rdot_3; rddot_3; rdot_rts; rddot_rts]; %a 24x1 matrix NomOrb (from 5.2.2. The Nominal Orbit) function [X_traj, VirtualOrb_Diff, Leg_Length, EE, E1, E2, E3] = NomOrb() %*********************************************************** %Converts Steve Hughes's set of orbital elements for an equilateral %formation to position and velocity % %REF.: Steve Hughes's paper "Preliminary Optimal ORbit Design for the %Laser %Interferometer Space Antenna (LISA)" % %CALLS: Orb2X2.m, X2Orb.m, LISA_EOMderivs.m %*********************************************************** %***CONSTANTS******************************************** mus = 1.32712428e11*((1e3)^3); %m3/s2, Sun's gravitational parameter d = 5e6*1000; %m, nominal leg length r2d = (180/pi); %radians to degrees converter Sun = [0 0 0]; %Sun is the Center %***ORBITAL ELEMENTS**************************** 165 ae = (1.49598023e8)*1000; %m, Semimajor Axis of Earth's orbit around Sun e = d/((2*sqrt(3))*ae); %unitless, Eccentricity i = (d/(2*ae))*r2d; %deg., Inclination (low inclination) w = 90; %90 deg. = (pi/2) or 270 deg. = (3*pi)/2 Argument of Periapsis Om1 = 0; %deg., Right Ascension of the Ascending Node of LISA 1 Om2 = Om1 + 120; %deg., 120 deg. = ((2*pi)/3) Right Asension of the Ascending Node of LISA 2 Om3 = Om1 - 120; %deg., Right Asension of the Ascending Node of LISA 3 M1 = 0; %deg., Mean Anomaly of LISA 1 M2 = M1 - 120; %deg., 120 deg = ((2*pi)/3) Mean Anomaly of LISA 2 M3 = M1 + 120; %deg., Mean Anomaly of LISA 3 %***CALL Orb2X2.m to convert to R & V************************ %Ephem4 = [a e i Om w M] %NOTE: Degrees will be converted to Radians in the code Orb2X2 EphemL1 = [ae e i Om1 w M1]; EphemL2 = [ae e i Om2 w M2]; EphemL3 = [ae e i Om3 w M3]; %[R, V] = Orb2X2(Ephem, t, mu, type) [R1, V1] = Orb2X2(EphemL1, 0, mus, 4); [R2, V2] = Orb2X2(EphemL2, 0, mus, 4); [R3, V3] = Orb2X2(EphemL3, 0, mus, 4); %***Position & Velocity of Earth Relative to Center of LISA*********** Rvs = [0 ae 0]; %1x3, 3-D initial Reference position of center of LISA's %virtual orbit around the Sun, lies all along the y-axis initially wvs = sqrt(mus/ae^3); %angular rate of virtual orbit v = wvs*ae; %m/s, tangential velocities of the spacecrafts Vvs = [-v 0 0]; %initial velocity of the virtual orbit in the x-axis th = 20*(pi/180); %Earth is 20 deg. rotation about the z-axis ahead of the %LISA formation converted to radians REt = [cos(th) -sin(th) 0; sin(th) cos(th) 0; 0 0 1]*Rvs'; %the transformed %position of the Earth VEt = [cos(th) -sin(th) 0; sin(th) cos(th) 0; 0 0 1]*Vvs'; %the transformed %velocity of the Earth %***CALL X2Orb.m to GET ORBITAL ELEMENTS OF LISA & %EARTH******************* %[E] = X2Orb(R, V, t, mu) %E = [a e i Om w M f Period] [EE] = X2Orb(REt, VEt, 0, mus); [E1] = X2Orb(R1, V1, 0, mus); [E2] = X2Orb(R2, V2, 0, mus); [E3] = X2Orb(R3, V3, 0, mus); %***LEG LENGTHS****************************************** L12 = (norm(R2 - R1)); %nx1, leg lengths L23 = (norm(R3 - R2)); L31 = (norm(R1 - R3)); Leg_Length = [L12, L23, L31]; %***DIFFERENCE FROM VIRTUAL ORBIT TO EACH LISA & %EARTH********************* L1C = norm(Rvs' - R1); L2C = norm(Rvs' - R2); L3C = norm(Rvs' - R3); EC2 = norm(Rvs' - REt); VirtualOrb_Diff = [L1C, L2C, L3C, EC2]; %***PROPAGATE THE ORBIT OVER THE YEAR USING THE %EOM************************ to = 0; I = (60*60)/2; %sec., 1/2 HOUR Intervals for the entire YEAR tf = 366.25*24*60*60; %sec, 1 SIDEREAL YEAR ts = 0:I:tf; 166 Xstar = [R1; V1; R2; V2; R3; V3; REt; VEt]; [t, Xs] = ode45('LISA_EOMderivs', ts, Xstar, odeset('RelTol',1e- 8,'AbsTol',1e-8)); X_traj = [ts(:) Xs(:,:)]; %n x 25 matrix LISA_STMderivs.m (from Section 5.2.3 The Algorithm) function stmdot = LISA_stmderive(t, Xstm) %*********************************************************** %function stmdot = LISA_stmderive(t, Xstm) % %Function that requires ode45 to integrate the state transition matrix over %a timespan % %INPUTS: t % Xstm = the state vector and the state transition matirx %OUTPUTS: stmdot % %Command Window Prompts %[t, Xstm] = ode45('LISA_stmderive', timespan, Xstmi, odeset('RelTol', %1e-8,'AbsTol',1e-8)) %*********************************************************** mus = (1.32712428e20); %m3/s2, Sun's gravitational parameter r1 = Xstm(1:3); %m, LISA spacecraft #1 location rmag_1 = norm(r1); rdot_1 = Xstm(4:6); %m/sec, LISA spacecraft #1 velocity rddot_1 = -(mus/(rmag_1^3))*r1; %m/s2, EOM for spacecraft #1 r2 = Xstm(7:9); %m, LISA spacecraft #2 location rmag_2 = norm(r2); rdot_2 = Xstm(10:12); %m/sec, LISA spacecraft #2 velocity rddot_2 = -(mus/(rmag_2^3))*r2; %m/s2, EOM for spacecraft #2 r3 = Xstm(13:15); %m, LISA spacecraft #3 location rmag_3 = norm(r3); rdot_3 = Xstm(16:18); %m/sec, LISA spacecraft #3 velocity rddot_3 = -(mus/(rmag_3^3))*r3; %m/s2, EOM for spacecraft #3 rts = Xstm(19:21); %m, tracking station location IS Earth (treat Earth as a point) rts_mag = norm(rts); rdot_rts = Xstm(22:24); %m/sec, Earth velocity rddot_rts = -(mus/(rts_mag^3))*rts; %m/s2, EOM for spacecraft #3 Xdot = [rdot_1; rddot_1; rdot_2; rddot_2; rdot_3; rddot_3; rdot_rts; rddot_rts]; %a 24x1 matrix %---A matrix, not constant because the variable r in the equations I = eye(3,3); Z3 = zeros(3,3); Gamma1 = -(mus/((norm(r1))^3))*I + 3*(((mus*r1)/(norm(r1)^5))*r1'); % 3x3, Gravity gradient matrix Gamma2= -(mus/((norm(r2))^3))*I + 3*(((mus*r2)/(norm(r2)^5))*r2'); % 3x3, Gravity gradient matrix Gamma3 = -(mus/((norm(r3))^3))*I + 3*(((mus*r3)/(norm(r3)^5))*r3'); % 3x3, Gravity gradient matrix Gammae = -(mus/((norm(rts))^3))*I + 3*(((mus*rts)/(norm(rts)^5))*rts'); % 3x3, Gravity gradient matrix A = [Z3 I Z3 Z3 Z3 Z3 Z3 Z3; Gamma1 Z3 Z3 Z3 Z3 Z3 Z3 Z3; Z3 Z3 Z3 I Z3 Z3 Z3 Z3;... Z3 Z3 Gamma2 Z3 Z3 Z3 Z3 Z3; Z3 Z3 Z3 Z3 Z3 I Z3 Z3; Z3 Z3 Z3 Z3 Gamma3 Z3 Z3 Z3;... Z3 Z3 Z3 Z3 Z3 Z3 Z3 I; Z3 Z3 Z3 Z3 Z3 Z3 Gammae Z3]; %24x24 matrix stm = reshape(Xstm(25:600), 24, 24); stmdot = A*stm; stmdot = [Xdot; (reshape(stmdot, 576, 1))]; %a 600x1 matrix 167 H_partials.m (from Section 5.2.3. The Algorithm) function [Hi_tilda, Hl_tilda, Hq_tilda] = H_partials7(Gi, GiR, Xstar, obs) %*********************************************************** %*********************************************************** X = Xstar(1,:); Z13 = zeros(1,3); %***THE PARTIALS %Range wrt Position dp1dx1 = (X(1)-X(19))/Gi(1); %m, L1 Range wrt L1 Position dp1dy1 = (X(2)-X(20))/Gi(1); dp1dz1 = (X(3)-X(21))/Gi(1); dp1dR = [dp1dx1 dp1dy1 dp1dz1]; dp2dx2 = (X(7)-X(19))/Gi(2); %m, L2 Range wrt L2 Position dp2dy2 = (X(8)-X(20))/Gi(2); dp2dz2 = (X(9)-X(21))/Gi(2); dp2dR = [dp2dx2 dp2dy2 dp2dz2]; dp3dx3 = (X(13)-X(19))/Gi(3); %m, L3 Range wrt L3 Position dp3dy3 = (X(14)-X(20))/Gi(3); dp3dz3 = (X(15)-X(21))/Gi(3); dp3dR = [dp3dx3 dp3dy3 dp3dz3]; %RANGE wrt VELOCITY dp1dV = Z13; dp2dV = Z13; dp3dV = Z13; %***RANGE RATE %LISA 1 RANGE RATE wrt its POSITION dpdot1dx1 = ((X(4)-X(22))/Gi(1)) - ((Gi(4)*(X(1)-X(19)))/(Gi(1)^2)); dpdot1dy1 = ((X(5)-X(23))/Gi(1)) - ((Gi(4)*(X(2)-X(20)))/(Gi(1)^2)); dpdot1dz1 = ((X(6)-X(24))/Gi(1)) - ((Gi(4)*(X(3)-X(21)))/(Gi(1)^2)); dpdot1dR = [dpdot1dx1 dpdot1dy1 dpdot1dz1]; %LISA 2 RANGE RATE wrt its POSITION dpdot2dx2 = ((X(10)-X(22))/Gi(2)) - ((Gi(5)*(X(7)-X(19)))/(Gi(2)^2)); dpdot2dy2 = ((X(11)-X(23))/Gi(2)) - ((Gi(5)*(X(8)-X(20)))/(Gi(2)^2)); dpdot2dz2 = ((X(12)-X(24))/Gi(2)) - ((Gi(5)*(X(9)-X(21)))/(Gi(2)^2)); dpdot2dR = [dpdot2dx2 dpdot2dy2 dpdot2dz2]; %LISA 3 RANGE RATE wrt its POSITION dpdot3dx3 = ((X(16)-X(22))/Gi(3)) - ((Gi(6)*(X(13)-X(19)))/(Gi(3)^2)); dpdot3dy3 = ((X(17)-X(23))/Gi(3)) - ((Gi(6)*(X(14)-X(20)))/(Gi(3)^2)); dpdot3dz3 = ((X(18)-X(24))/Gi(3)) - ((Gi(6)*(X(15)-X(21)))/(Gi(3)^2)); dpdot3dR = [dpdot3dx3 dpdot3dy3 dpdot3dz3]; %LISA 1,2,3 Range Rate wrt its VELOCITY dpdot1dV = [dp1dx1 dp1dy1 dp1dz1]; dpdot2dV = [dp2dx2 dp2dy2 dp2dz2]; dpdot3dV = [dp3dx3 dp3dy3 dp3dz3]; %***DECLINATION PARTIALS %LISA 1 declination wrt its POSITION dd1dx1 =((X(3)-X(21))* (-X(1) + X(19)))/... ((Gi(1)^3)*(1-(((X(3)-X(21))^2)/(Gi(1)^2)))^(1/2)); dd1dy1 = ((X(3)-X(21))* (-X(2) + X(20)))/... ((Gi(1)^3)*(1-(((X(3)-X(21))^2)/(Gi(1)^2)))^(1/2)); 168 dd1dz1 = ((1/Gi(1)) - (.5*(((X(3)-X(21))*((2*X(3)) - (2*X(21))))/(Gi(1)^3))))/... (1-(((X(3)-X(21))^2)/(Gi(1)^2)))^(1/2); dd1dR = [dd1dx1 dd1dy1 dd1dz1]; %LISA 2 declination wrt its POSITION dd2dx2 = ((X(9)-X(21))* (-X(7) + X(19)))/... ((Gi(2)^3)*(1-(((X(9)-X(21))^2)/(Gi(2)^2)))^(1/2)); dd2dy2 = ((X(9)-X(21))* (-X(8) + X(20)))/... ((Gi(2)^3)*(1-(((X(9)-X(21))^2)/(Gi(2)^2)))^(1/2)); dd2dz2 = ((1/Gi(2)) - (.5*(((X(9)-X(21))*((2*X(9)) - (2*X(21))))/(Gi(2)^3))))/... (1-(((X(9)-X(21))^2)/(Gi(2)^2)))^(1/2); dd2dR = [dd2dx2 dd2dy2 dd2dz2]; %LISA 3 declination wrt its POSITION dd3dx3 = ((X(15)-X(21))* (-X(13) + X(19)))/... ((Gi(3)^3)*(1-(((X(15)-X(21))^2)/(Gi(3)^2)))^(1/2)); dd3dy3 = ((X(15)-X(21))* (-X(14) + X(20)))/... ((Gi(3)^3)*(1-(((X(15)-X(21))^2)/(Gi(3)^2)))^(1/2)); dd3dz3 = ((1/Gi(3)) - (.5*(((X(15)-X(21))*((2*X(15)) - (2*X(21))))/(Gi(3)^3))))/... (1-(((X(15)-X(21))^2)/(Gi(3)^2)))^(1/2); dd3dR = [dd3dx3 dd3dy3 dd3dz3]; %***Delta wrt VELOCITY dd1dV = Z13; dd2dV = Z13; dd3dV = Z13; %RIGHT ASCENSION PARTIALS %LISA 1 right ascension wrt its position da1dx1 = (-1*(X(2)-X(20)))/(((X(1)-X(19))^2)*(1 + ((X(2)- X(20))^2)/((X(1)-X(19))^2))); da1dy1 = 1/((X(1)-X(19))*(1 + ((X(2)-X(20))^2)/((X(1)-X(19))^2))); da1dz1 = 0; da1dR = [da1dx1 da1dy1 da1dz1]; %LISA 2 right ascension wrt its position da2dx2 = (-1*(X(8)-X(20)))/(((X(7)-X(19))^2)*(1+ ((X(8)- X(20))^2)/((X(7)-X(19))^2))); da2dy2 = 1/((X(7)-X(19))*(1 + ((X(8)-X(20))^2)/((X(7)-X(19))^2))); da2dz2 = 0; da2dR = [da2dx2 da2dy2 da2dz2]; %LISA 3 right ascension wrt its position da3dx3 = (-1*(X(14)-X(20)))/(((X(13)-X(19))^2)*(1+ ((X(14)- X(20))^2)/((X(13)-X(19))^2))); da3dy3 = 1/((X(13)-X(19))*(1 + ((X(14)-X(20))^2)/((X(13)-X(19))^2))); da3dz3 = 0; da3dR = [da3dx3 da3dy3 da3dz3]; %***Alpha wrt VELOCITY da1dV = Z13; da2dV = Z13; da3dV = Z13; %***INTERSPACECRAFT RANGE AND RANGE RATE %*********************************************************** %LISA 1 range wrt LISA 3 position dpL13dx1 = (X(1) - X(13))/GiR(1); dpL13dy1 = (X(2) - X(14))/GiR(1); dpL13dz1 = (X(3) - X(15))/GiR(1); dpL13dR1 = [dpL13dx1 dpL13dy1 dpL13dz1]; dpL13dV1 = Z13; %LISA 1 range rate wrt LISA 3 position 169 dpdotL13dx1 = ((X(4)-X(16))/GiR(1)) - ((GiR(4)*(X(1)- X(13)))/(GiR(1)^2)); dpdotL13dy1 = ((X(5)-X(17))/GiR(1)) - ((GiR(4)*(X(2)- X(14)))/(GiR(1)^2)); dpdotL13dz1 = ((X(6)-X(18))/GiR(1)) - ((GiR(4)*(X(3)- X(15)))/(GiR(1)^2)); dpdotL13dR1 =[dpdotL13dx1 dpdotL13dy1 dpdotL13dz1]; dpdotL13dV1 = [dpL13dx1 dpL13dy1 dpL13dz1]; %LISA 3 range wrt LISA 1 position dpL13dx3 = (-X(1) + X(13))/GiR(1); dpL13dy3 = (-X(2) + X(14))/GiR(1); dpL13dz3 = (-X(3) + X(15))/GiR(1); dpL13dR3 = [dpL13dx3 dpL13dy3 dpL13dz3]; dpL13dV3 = Z13; %LISA 3 range rate wrt LISA 1 position dpdotL13dx3 = ((-X(4)+X(16))/GiR(1)) - ((GiR(4)*(- X(1)+X(13)))/(GiR(1)^2)); dpdotL13dy3 = ((-X(5)+X(17))/GiR(1)) - ((GiR(4)*(- X(2)+X(14)))/(GiR(1)^2)); dpdotL13dz3 = ((-X(6)+X(18))/GiR(1)) - ((GiR(4)*(- X(3)+X(15)))/(GiR(1)^2)); dpdotL13dR3 =[dpdotL13dx3 dpdotL13dy3 dpdotL13dz3]; dpdotL13dV3 = [dpL13dx3 dpL13dy3 dpL13dz3]; %*********************************************************** %LISA 2 wrt LISA 1 %LISA 2 range wrt LISA 1 position dpL21dx2 = (X(7) - X(1))/GiR(2); dpL21dy2 = (X(8) - X(2))/GiR(2); dpL21dz2 = (X(9) - X(3))/GiR(2); dpL21dR2 = [dpL21dx2 dpL21dy2 dpL21dz2]; dpL21dV2 = Z13; %LISA 2 range rate wrt LISA 1 position dpdotL21dx2 = ((X(10)-X(4))/GiR(2)) - ((GiR(5)*(X(7)- X(1)))/(GiR(2)^2)); dpdotL21dy2 = ((X(11)-X(5))/GiR(2)) - ((GiR(5)*(X(8)- X(2)))/(GiR(2)^2)); dpdotL21dz2 = ((X(12)-X(6))/GiR(2)) - ((GiR(5)*(X(9)- X(3)))/(GiR(2)^2)); dpdotL21dR2 =[dpdotL21dx2 dpdotL21dy2 dpdotL21dz2]; dpdotL21dV2 = [dpL21dx2 dpL21dy2 dpL21dz2]; %LISA 1 range wrt LISA 2 position dpL21dx1 = (-X(7) + X(1))/GiR(2); dpL21dy1 = (-X(8) + X(2))/GiR(2); dpL21dz1 = (-X(9) + X(3))/GiR(2); dpL21dR1 = [dpL21dx1 dpL21dy1 dpL21dz1]; dpL21dV1 = Z13; %LISA 1 range rate wrt LISA 2 position dpdotL21dx1 = ((-X(10)+X(4))/GiR(2)) - ((GiR(5)*(- X(7)+X(1)))/(GiR(2)^2)); dpdotL21dy1 = ((-X(11)+X(5))/GiR(2)) - ((GiR(5)*(- X(8)+X(2)))/(GiR(2)^2)); dpdotL21dz1 = ((-X(12)+X(6))/GiR(2)) - ((GiR(5)*(- X(9)+X(3)))/(GiR(2)^2)); dpdotL21dR1 =[dpdotL21dx1 dpdotL21dy1 dpdotL21dz1]; dpdotL21dV1 = [dpL21dx1 dpL21dy1 dpL21dz1]; %*********************************************************** 170 %LISA 3 wrt LISA 2 %LISA 3 range wrt LISA 2 position dpL32dx3 = (X(13) - X(7))/GiR(3); dpL32dy3 = (X(14) - X(8))/GiR(3); dpL32dz3 = (X(15) - X(9))/GiR(3); dpL32dR3 = [dpL32dx3 dpL32dy3 dpL32dz3]; dpL32dV3 = Z13; %LISA 3 range rate wrt LISA 2 position dpdotL32dx3 = ((X(16)-X(10))/GiR(3)) - ((GiR(6)*(X(13)- X(7)))/(GiR(3)^2)); dpdotL32dy3 = ((X(17)-X(11))/GiR(3)) - ((GiR(6)*(X(14)- X(8)))/(GiR(3)^2)); dpdotL32dz3 = ((X(18)-X(12))/GiR(3)) - ((GiR(6)*(X(15)- X(9)))/(GiR(3)^2)); dpdotL32dR3 =[dpdotL32dx3 dpdotL32dy3 dpdotL32dz3]; dpdotL32dV3 = [dpL32dx3 dpL32dy3 dpL32dz3]; %LISA 2 range wrt LISA 3 position dpL32dx2 = (-X(13) + X(7))/GiR(3); dpL32dy2 = (-X(14) + X(8))/GiR(3); dpL32dz2 = (-X(15) + X(9))/GiR(3); dpL32dR2 = [dpL32dx2 dpL32dy2 dpL32dz2]; dpL32dV2 = Z13; %LISA 2 range rate wrt LISA 3 position dpdotL32dx2 = ((-X(16)+X(10))/GiR(3)) - ((GiR(6)*(- X(13)+X(7)))/(GiR(3)^2)); dpdotL32dy2 = ((-X(17)+X(11))/GiR(3)) - ((GiR(6)*(- X(14)+X(8)))/(GiR(3)^2)); dpdotL32dz2 = ((-X(18)+X(12))/GiR(3)) - ((GiR(6)*(- X(15)+X(9)))/(GiR(3)^2)); dpdotL32dR2 =[dpdotL32dx2 dpdotL32dy2 dpdotL32dz2]; dpdotL32dV2 = [dpL32dx2 dpL32dy2 dpL32dz2]; switch(obs) case 0 % DSN & VLBI %Hi, Hl, Hq into 4x18 matrices [1x3, 1x3, 1x3, 1x3, 1x3, 1x3] Hi_tilda = [dp1dR dp1dV Z13 Z13 Z13 Z13;... dpdot1dR dpdot1dV Z13 Z13 Z13 Z13;... dd1dR dd1dV Z13 Z13 Z13 Z13;... da1dR da1dV Z13 Z13 Z13 Z13]; Hl_tilda = [Z13 Z13 dp2dR dp2dV Z13 Z13;... Z13 Z13 dpdot2dR dpdot2dV Z13 Z13;... Z13 Z13 dd2dR dd2dV Z13 Z13;... Z13 Z13 da2dR da2dV Z13 Z13]; Hq_tilda = [Z13 Z13 Z13 Z13 dp3dR dp3dV;... Z13 Z13 Z13 Z13 dpdot3dR dpdot3dV;... Z13 Z13 Z13 Z13 dd3dR dd3dV;... Z13 Z13 Z13 Z13 da3dR da3dV]; case 1 %DSN, VLBI, RELATIVE %H 18x18 matrix [1x3, 1x3, 1x3, 1x3, 1x3, 1x3] Hi_tilda = [dp1dR dp1dV Z13 Z13 Z13 Z13;... dpdot1dR dpdot1dV Z13 Z13 Z13 Z13;... dd1dR dd1dV Z13 Z13 Z13 Z13;... da1dR da1dV Z13 Z13 Z13 Z13;... dpL13dR1 dpL13dV1 Z13 Z13 dpL13dR3 dpL13dV3;... dpdotL13dR1 dpdotL13dV1 Z13 Z13 dpdotL13dR3 dpdotL13dV3;...%wrt LISA 1 Z13 Z13 dp2dR dp2dV Z13 Z13;... Z13 Z13 dpdot2dR dpdot2dV Z13 Z13;... Z13 Z13 dd2dR dd2dV Z13 Z13;... Z13 Z13 da2dR da2dV Z13 Z13;... dpL21dR1 dpL21dV1 dpL21dR2 dpL21dV2 Z13 Z13;... 171 dpdotL21dR1 dpdotL21dV1 dpdotL21dR2 dpdotL21dV2 Z13 Z13;...%wrt LISA 2 Z13 Z13 Z13 Z13 dp3dR dp3dV;... Z13 Z13 Z13 Z13 dpdot3dR dpdot3dV;... Z13 Z13 Z13 Z13 dd3dR dd3dV;... Z13 Z13 Z13 Z13 da3dR da3dV;... Z13 Z13 dpL32dR2 dpL32dV2 dpL32dR3 dpL32dV3;... Z13 Z13 dpdotL32dR2 dpdotL32dV2 dpdotL32dR3 dpdotL32dV3]; %wrt LISA 3 Hl_tilda = [Z13]; Hq_tilda = [Z13]; case 2 %DSN only %Hi, Hl, Hq into 2x18 matrices [1x3, 1x3, 1x3, 1x3, 1x3, 1x3] Hi_tilda = [dp1dR dp1dV Z13 Z13 Z13 Z13;... dpdot1dR dpdot1dV Z13 Z13 Z13 Z13]; Hl_tilda = [Z13 Z13 dp2dR dp2dV Z13 Z13;... Z13 Z13 dpdot2dR dpdot2dV Z13 Z13]; Hq_tilda = [Z13 Z13 Z13 Z13 dp3dR dp3dV;... Z13 Z13 Z13 Z13 dpdot3dR dpdot3dV]; case 3 %DSN & RELATIVE %H 12x18 matrix [1x3, 1x3, 1x3, 1x3, 1x3, 1x3] Hi_tilda = [dp1dR dp1dV Z13 Z13 Z13 Z13;... dpdot1dR dpdot1dV Z13 Z13 Z13 Z13;... dpL13dR1 dpL13dV1 Z13 Z13 dpL13dR3 dpL13dV3;... dpdotL13dR1 dpdotL13dV1 Z13 Z13 dpdotL13dR3 dpdotL13dV3;...%wrt LISA 1 Z13 Z13 dp2dR dp2dV Z13 Z13;... Z13 Z13 dpdot2dR dpdot2dV Z13 Z13;... dpL21dR1 dpL21dV1 dpL21dR2 dpL21dV2 Z13 Z13;... dpdotL21dR1 dpdotL21dV1 dpdotL21dR2 dpdotL21dV2 Z13 Z13;...%wrt LISA 2 Z13 Z13 Z13 Z13 dp3dR dp3dV;... Z13 Z13 Z13 Z13 dpdot3dR dpdot3dV;... Z13 Z13 dpL32dR2 dpL32dV2 dpL32dR3 dpL32dV3;... Z13 Z13 dpdotL32dR2 dpdotL32dV2 dpdotL32dR3 dpdotL32dV3]; %wrt LISA 3 Hl_tilda = [Z13]; Hq_tilda = [Z13]; otherwise error('Invalid case (0 <= obs <= 3)'); end RealObs.m (from Section 5.2.4 Simulation of Measurement Noise) function [Gi, L1tot, L2tot, L3tot] = RealObs(to, tf, I, X, Xref, sigma_range, sigma_rrate, sigma_alpha, sigma_delta, sigma_ci, sigma_f, B, plt) %*********************************************************** %This function calculates the "perfect" observations and also calculates %the corresponding "noisy" measurements. % %INPUT : to = Initial time to start propagation (sec) % tf = Final time to end propagation (sec) % I = Time intervals (sec) % XO = An Initial Condition of the State (1x24) % sigma_range (m), sigma_rrate(m/s), sigma_alpha(rad), sigma_delta(rad) = The noises % associated with each measurement % br(m), brr(m/s), ba(rad), bs(rad) = biases in measurements % plt = plot command =0 "off", = 1 "on" % %OUTPUT : Gi = 1x12 vector of calculated range, range rate, alpha and delta 172 % angles of each LISA % L1tot, L2tot, L3tot = 4x1 vector of Measurement Deviations [range, range rate, alpha and delta % angles] of each LISA %*********************************************************** %---INITIALIZATION if nargin < 13, plt = 0; %default to turn plot command off end %***Parse Truth vector (X)into position & velocity components************** R1 = X(1:3); %m (1x3) R2 = X(7:9); R3 = X(13:15); R = [R1 R2 R3]; %(1x9) V1 = X(4:6); %m/s (1x3) V2 = X(10:12); V3 = X(16:18); V = [V1 V2 V3]; %(1x9) Rts= X(19:21); %(1x3) Vts = X(22:24); %***Parse Reference vector (Xref) into position & velocity components****** R1s = Xref(1:3); %m (1x3) R2s = Xref(7:9); R3s = Xref(13:15); Rs = [R1s R2s R3s]; %(1x9) V1s = Xref(4:6); %m/s (1x3) V2s = Xref(10:12); V3s = Xref(16:18); Vs = [V1s V2s V3s]; %(1x9) Rts_ref = Xref(19:21); %(1x3) Vts_ref = Xref(22:24); %***CALCULATE THE "ACTUAL" MEASURMENTS (i.e. Truth plus noise and biases)** %RANGE p1 = norm(R1 - Rts); %m (1x1) range scalar for LISA 1 p2 = norm(R2 - Rts); p3 = norm(R3 - Rts); %RANGE RATE pdot1 = (((R1(1) - Rts(1))*(V1(1) - Vts(1))) + ((R1(2) - Rts(2))*... (V1(2) - Vts(2))) + ((R1(3) - Rts(3))*(V1(3) - Vts(3))))/p1; %m/s pdot2 =(((R2(1) - Rts(1))*(V2(1)- Vts(1))) + ((R2(2) - Rts(2))*... (V2(2) - Vts(2))) + ((R2(3) - Rts(3))*(V2(3) - Vts(3))))/p2; pdot3 =(((R3(1) - Rts(1))*(V3(1)- Vts(1))) + ((R3(2) - Rts(2))*... (V3(2) - Vts(2))) + ((R3(3) - Rts(3))*(V3(3) - Vts(3))))/p3; %ALPHA alpha1 = atan2((R1(2)- Rts(2)),(R1(1)-Rts(1))); %rad (1x1) alpha2 = atan2((R2(2)- Rts(2)),(R2(1)-Rts(1))); alpha3 = atan2((R3(2)- Rts(2)),(R3(1)-Rts(1))); %DELTA delta1 = asin((R1(3)- Rts(3))/p1); %rad (1x1) delta2 = asin((R2(3)- Rts(3))/p2); delta3 = asin((R3(3)- Rts(3))/p3); %***ADD NOISE TO THE TRUTH************************************************* Yi_1 = p1 + sigma_range(1) + sigma_ci(1) + B(1); %m Yl_1 = p2 + sigma_range(2) + sigma_ci(2) + B(2); Yq_1 = p3 + sigma_range(3) + sigma_ci(3) + B(3); Yi_2 = (pdot1) + sigma_rrate(1) + sigma_f(1); %m/s 173 Yl_2 = (pdot2) + sigma_rrate(2) + sigma_f(2); Yq_2 = (pdot3) + sigma_rrate(3) + sigma_f(3); Yi_3 = (alpha1) + sigma_alpha(1); %rad Yl_3 = (alpha2) + sigma_alpha(2); Yq_3 = (alpha3) + sigma_alpha(3); Yi_4 = (delta1) + sigma_delta(1); %rad Yl_4 = (delta2) + sigma_delta(2); Yq_4 = (delta3) + sigma_delta(3); %***CALCULATE THE MEASURMENTS (i.e. Reference Observations)**************** %RANGE + ADD KNOWN BIASES TO CALCULATED MEASUREMENT p1s = norm(R1s - Rts_ref);% + B(1); %m (1x1) range scalar for LISA 1 p2s = norm(R2s - Rts_ref);% + B(2); p3s = norm(R3s - Rts_ref);% + B(3); %RANGE RATE pdot1s = ((((R1s(1) - Rts_ref(1))*(V1s(1) - Vts_ref(1))) + ((R1s(2) - Rts_ref(2))*... (V1s(2) - Vts_ref(2))) + ((R1s(3) - Rts_ref(3))*(V1s(3) - Vts_ref(3))))/p1s); %m/s pdot2s =((((R2s(1) - Rts_ref(1))*(V2s(1)- Vts_ref(1))) + ((R2s(2) - Rts_ref(2))*... (V2s(2) - Vts_ref(2))) + ((R2s(3) - Rts_ref(3))*(V2s(3) - Vts_ref(3))))/p2s); pdot3s =((((R3s(1) - Rts_ref(1))*(V3s(1)- Vts_ref(1))) + ((R3s(2) - Rts_ref(2))*... (V3s(2) - Vts_ref(2))) + ((R3s(3) - Rts_ref(3))*(V3s(3) - Vts_ref(3))))/p3s); %ALPHA alpha1s = (atan2((R1s(2) - Rts_ref(2)),(R1s(1) - Rts_ref(1)))); %rad (1x1) alpha2s = (atan2((R2s(2) - Rts_ref(2)),(R2s(1) - Rts_ref(1)))); alpha3s = (atan2((R3s(2) - Rts_ref(2)),(R3s(1) - Rts_ref(1)))); %DELTA delta1s = (asin((R1s(3) - Rts_ref(3))/p1s)); %rad (1x1) delta2s = (asin((R2s(3) - Rts_ref(3))/p2s)); delta3s = (asin((R3s(3) - Rts_ref(3))/p3s)); %*********************************************************** Yi = [Yi_1' Yl_1' Yq_1' Yi_2' Yl_2' Yq_2' Yi_3' Yl_3' Yq_3' Yi_4' Yl_4' Yq_4']; %Real Measurements (1x12) Gi = [p1s' p2s' p3s' pdot1s' pdot2s' pdot3s' alpha1s' alpha2s' alpha3s' delta1s' delta2s' delta3s']; %Calculated Measurements (1x12) diff = Yi-Gi; %---COMPUTE MEASUREMENT DEVIATIONS, yi (1x1), AT TIME ti %LISA 1, 2 & 3 measurement deviations; measurement(Yi) - Reference %Measurement(Gi) yi_1 = Yi(1) - Gi(1); % (1x1), range measurement deviations yl_1 = Yi(2) - Gi(2); yq_1 = Yi(3) - Gi(3); yi_2 = Yi(4) - Gi(4); %range rate measurement deviations yl_2 = Yi(5) - Gi(5); yq_2 = Yi(6) - Gi(6); yi_3 = Yi(10) - Gi(10); %delta measurement deviations yl_3 = Yi(11) - Gi(11); yq_3 = Yi(12) - Gi(12); yi_4 = Yi(7) - Gi(7); %alpha measurement deviations yl_4 = Yi(8) - Gi(8); yq_4 = Yi(9) - Gi(9); 174 %Total Measurement deviations L1tot = [yi_1; yi_2; yi_3; yi_4]; %(4x1) L2tot = [yl_1; yl_2; yl_3; yl_4]; %(4x1) L3tot = [yq_1; yq_2; yq_3; yq_4]; %(4x1) %PLOTS (x,y,z) if plt == 1, figure subplot(2,2,1) plot(timespan, Gi(:,1),'b-', timespan, Gi(:,2),'g--', timespan, Gi(:,3), 'r-.',... timespan, Yi(:,1),'b*', timespan, Yi(:,2),'g>', timespan, Yi(:,3),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Range(m)') legend('LISA 1 (calc.)', 'LISA 2(calc.)', 'LISA 3(calc.)','LISA 1(real)', 'LISA 2(real)', 'LISA 3(real)') title('LISA Range vs. Time') grid on subplot(2,2,2) plot(timespan, Gi(:,4),'b-', timespan, Gi(:,5), 'g--', timespan, Gi(:,6),'r-.',... timespan, Yi(:,4),'b*', timespan, Yi(:,5),'b>', timespan, Yi(:,6),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Range Rate(m/sec)') %legend('LISA 1', 'LISA 2', 'LISA 3') title('LISA Range Rate vs. Time') grid on subplot(2,2,3) plot(timespan, Gi(:,7),'b-', timespan, Gi(:,8), 'g--', timespan, Gi(:,9),'r-.',... timespan, Yi(:,7),'b*', timespan, Yi(:,8),'g>', timespan, Yi(:,9),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Alpha angle (rad)') %legend('LISA 1', 'LISA 2', 'LISA 3') title('LISA Alpha Angle vs. Time') grid on subplot(2,2,4) plot(timespan, Gi(:,10),'b-', timespan, Gi(:,11), 'g--', timespan, Gi(:,12),'r- .',... timespan, Yi(:,10),'b*', timespan, Yi(:,11),'g>', timespan, Yi(:,12),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Delta angle (rad)') %legend('LISA 1', 'LISA 2', 'LISA 3') title('LISA Delta Angle vs. Time') grid on else end InterSCD.m (from Section 5.2.4 Simulation of Measurement Noise) function [GiR, L13tot, L21tot, L32tot] = InterSCD(to, tf, I, X, Xref, sigma_irange, sigma_irate, sigma_ciL, B, plt) %*********************************************************** %Program to generate the interspacecraft data for LISA. The range and range %rate information of each spacecraft wrt each other % %INPUT : to, tf, I = initial time, final time and Interval (sec) % XO = initial condition for state (1x24) % sigma_irange, sigma_irate = noises associated with the interspacecraft range and range rate % bir, birr = biases associated with the interspacecraft range and % range rate % plt = plot command, if 0 then "no plots", of 1 then "plot" % %OUTPUT : T_ISC = [pL13 pL21 pL32 pdotL13 pdotL21 pdotL32] % N_ISC = [Yi_5 Yl_5 Yq_5 Yi_6 Yl_6 Yq_6] % pL13 = range of LISA 1 wrt 3, with noise(and biases) = Yi_5 % pL21 = range of LISA 2 wrt 1 = Yl_5 % pL32 = range of LISA 3 wrt 2 = Yq_5 % pdotL13 = range rate of LISA 1 wrt 3 = Yi_6 175 % pdotL21 = range rate of LISA 2 wrt 1 = Yl_6 % pdotL32 = range rate of LISA 3 wrt 2 = Yq_6 % %*********************************************************** %---INITIALIZATION if nargin < 10, plt = 0; %default to turn plot command off end %timespan =(to:I:tf)'; %with I sec intervals of time %***PARSE TRUTH VECTOR(X) into Position & Velocity*********** R1 = X(1:3); %m (1x3) V1 = X(4:6); %m/s (1x3) R2 = X(7:9); V2 = X(10:12); R3 = X(13:15); V3 = X(16:18); RE = X(19:21); VE = X(22:24); %***PARSE REFERENCE VECTOR(Xref) into Position & Velocity R1s = Xref(1:3); %m (1x3) V1s = Xref(4:6); %m/s (1x3) R2s = Xref(7:9); V2s = Xref(10:12); R3s = Xref(13:15); V3s = Xref(16:18); REs = Xref(19:21); VEs = Xref(22:24); %***CALCULATE THE ACTUAL MEASUREMENTS(I.e. Truth plus %noise and biases)**** %RANGE pL13 = norm(R1-R3); %m (1x1) pL13_vec =(R1-R3); pL21 = norm(R2-R1); pL21_vec = (R2-R1); pL32 = norm(R3-R2); pL32_vec =(R3-R2); %RANGE RATE pdotL13 = (((R1(1)-R3(1))*(V1(1)- V3(1))) + ((R1(2)-R3(2))*... (V1(2) - V3(2))) + ((R1(3) - R3(3))*(V1(3) - V3(3))))/pL13; %m/s pdotL21 = (((R2(1)-R1(1))*(V2(1)- V1(1))) + ((R2(2)-R1(2))*... (V2(2) - V1(2))) + ((R2(3) - R1(3))*(V2(3)-V1(3))))/pL21; pdotL32 = (((R3(1)-R2(1))*(V3(1)- V2(1))) + ((R3(2)-R2(2))*... (V3(2) - V2(2))) + ((R3(3) - R2(3))*(V3(3)-V2(3))))/pL32; %***ADD THE NOISE TO THE TRUTH******************* % [GiR, YiR] = InterSCD2(to, ts(i), I, X, Xstar, S(i,13:15), S(i,16:18), %B(4:6), 0); %interspacecraft range stays in m Yi_5 = pL13 + sigma_irange(1) + sigma_ciL(1) + B(1); %m Yl_5 = pL21 + sigma_irange(2) + sigma_ciL(2) + B(2); Yq_5 = pL32 + sigma_irange(3) + sigma_ciL(3) + B(3); %keep in m/s Yi_6 = (pdotL13) + sigma_irate(1); %mm/s Yl_6 = (pdotL21) + sigma_irate(2); Yq_6 = (pdotL32) + sigma_irate(3); %***CALCULATE THE MEASUREMENTS i.e. Reference Observations %RANGE PLUS ADD KNOWN BIASES TO CALCULATED %MEASUREMENT pL13s = norm(R1s-R3s);%+ B(1); %m; %m (1x1) pL21s = norm(R2s-R1s);%+ B(2); %m; 176 pL32s = norm(R3s-R2s);%+ B(3); %m; %RANGE RATE pdotL13s = ((((R1s(1)-R3s(1))*(V1s(1)- V3s(1))) + ((R1s(2)-R3s(2))*... (V1s(2) - V3s(2))) + ((R1s(3) - R3s(3))*(V1s(3) - V3s(3))))/pL13s); %m/s pdotL21s = ((((R2s(1)-R1s(1))*(V2s(1)- V1s(1))) + ((R2s(2)-R1s(2))*... (V2s(2) - V1s(2))) + ((R2s(3) - R1s(3))*(V2s(3)-V1s(3))))/pL21s); pdotL32s = ((((R3s(1)-R2s(1))*(V3s(1)- V2s(1))) + ((R3s(2)-R2s(2))*... (V3s(2) - V2s(2))) + ((R3s(3) - R2s(3))*(V3s(3)-V2s(3))))/pL32s); %*********************************************************** YiR = [Yi_5' Yl_5' Yq_5' Yi_6' Yl_6' Yq_6']; %Real Measurements (1x6) GiR = [pL13s' pL21s' pL32s' pdotL13s' pdotL21s' pdotL32s']; %Calculated Measurements (1x6) %---COMPUTE MEASUREMENT DEVIATIONS, yi (1x1), AT TIME ti %LISA 1, 2 & 3 measurement deviations; measurement(Yi) ? Reference Measurement(Gi) yi_1 = YiR(1) - GiR(1); %(1x1), range measurement deviations yl_1 = YiR(2) - GiR(2); yq_1 = YiR(3) - GiR(3); yi_2 = YiR(4) - GiR(4); %range rate measurement deviations yl_2 = YiR(5) - GiR(5); yq_2 = YiR(6) - GiR(6); %Total Measurement deviations L13tot = [yi_1; yi_2]; %(2x1) L21tot = [yl_1; yl_2];%(2x1) L32tot = [yq_1; yq_2]; %(2x1) %PLOTS (x,y,z) if plt == 1, figure subplot(2,1,1) plot(timespan, GiR(:,1),'b-', timespan, GiR(:,2),'g--', timespan, GiR(:,3), 'r- .',... timespan, YiR(:,1),'b*', timespan, YiR(:,2),'g>', timespan, YiR(:,3),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Range(m)') legend('Leg 13', 'Leg 21', 'Leg 32','Leg 13', 'Leg 21', 'Leg 32') title('LISA Interspacecraft Range vs. Time') grid on subplot(2,1,2) plot(timespan, GiR(:,4),'b-', timespan, GiR(:,5), 'g--', timespan, GiR(:,6),'r- .',... timespan, YiR(:,4),'b*', timespan, YiR(:,5),'g>', timespan, YiR(:,6),'rx'); xlabel('Time (sec)') ylabel('Calculated / Real(noisy) Range Rate(m/sec)') title('LISA Interspacecraft Range Rate vs. Time') grid on else end LISA_WBLS.m (from Section 5.4 Orbit Accuracy) function [Xest, ABS, REL, RMS, var, PU_L1, PU_L2, PU_L3, err, erv, niter, Pcov] =MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R1, R2, R3, S, B, plt, obs) %*********************************************************** ******* %[Xest, DXr, DXv, DLr, DLv] = LISA_bp(to, tf, I, R1, R2, R3, AP, Po, xo_bar, %plt,obs) % %This function is a general Batch Processsor that takes in the %absolute spacecraft data (range, range rate, right ascension and declination) %and relative spacecraft data (interspacecraft range and range rate) %It can also be adapted to take in a priori information and weights. Error %analysis will also be done and graphed. 177 % %INPUTS : to_0 = initial starting time of data (sec) % tf = final time of data (sec) % I = time interval to space measurements (sec) % Xrefstart_0 = The Initial Guess Vector, or Reference % Start Vector (1x18), may be off by some errorstart from the % Xo_0 % Xo_0 = The Nominal Vector (aka Truth) (1x18) % R1, R2, R3 = each are 1xn array of weights associated with % each observation type for each LISA. Depends on size of % observation case % S = sigmas associated with each measurement % B = read in biases % plt= 0, "no" plot; =1, "yes" plot command is on % %Observation Cases: % obs = 0, DSN & VLBI % = 1, DSN, VLBI & RELATIVE % = 2, DSN only % = 3, DSN & Relative % %OUTPUTS: Xest = The final Batch Processor estimate of the state % ABS = The Absolute Error in position & velocity of the LISA % constellation % REL = The Relative Error in position & velocity of the LISA % constellation % RMS = The Root-Mean-Square of each LISA % var = The Sample Variance of each LISA % PU_L1, PU_L2, PU_L3 = Parameter Uncertainites % err = Position Error in each LISA % erv = Velocity Error in each LISA % %CALLS: LISA_EOMderivs.m, LISA_stmderive.m, RealObs.m, InterSCD.m, %H_partials.m % %REFERENCE: Fig. 4.6.1 (pp. 196-197) from B.Tapley's "Statistical Orbit %Determination" %*********************************************************** ******* %---0. INITIALIZATION format long e if nargin < 12 obs = 0; %default, takes into account all absolute (DSN & VLBI) observation types end if nargin <11 plt = 0; %default no plots end %---SET UP TOLERANCE FOR CONVERGENCE switch(obs) case 0, %DSN & VLBI tol = .03; case 1, %DSN, VLBI & RELATIVE tol = .03; case 2, %DSN tol = .3; %.3; %.2; %case 2 won't converge with lower like .09 case 3, %DSN & RELATIVE tol = .03; otherwise, error('Invalid case (0 <= obs <= 3)'); end RMSold = [1 1 1]; dRMS = 1; niter = 0; %set up number of iterations to = to_0; %initial time ts = to:I:tf; %(1xn) set up timespan Xo = Xo_0; 178 Xtru = Xo; %Nominal (or Truth) Vector Xk = Xrefstart_0; %Initial Guess Vector (i.e. Reference Vector) Xrefstart = Xk; %Will get updated from initial with each iteration in the program %---A. OUTER LOOP while dRMS >= tol %---INITIAL CONDITION FOR THE STATE Xstaro = Xrefstart; % (1x24) I.C. for integration = Xstar(ti-1) to = to_0; Xo = Xo_0; diff = Xo - Xo_0; %---DEFINE & GIVE THE INITIAL CONDITION FOR THE STATE %TRANSITION MATRIX %(STM) INTEGRATOR stmi = eye(24); %(24x24), identity matrix rstmi = reshape(stmi,576,1); %(576x1) Xstmi = [Xstaro'; rstmi]; %(600x1) I.C. for STM integration %---DEFINE THE M = Information Matrix, inverse of the covariance % & N = Reduced Residual Vector %Both matrices are different depending on if a priori knowledge switch(obs) case 0, %DSN & VLBI M1 = zeros(6,6); % (6x6) position & velocity of LISA 1 N1 = zeros(6,1); % (6x1) M2 = zeros(6,6); %LISA 2 N2 = zeros(6,1); M3 = zeros(6,6); %LISA 3 N3 = zeros(6,1); case 1, %DSN, VLBI & RELATIVE M = zeros(18,18); %LISA 1, 2 & 3 Postion & Velocity (18x18) N = zeros(18,1); %(18x1) case 2, %DSN M1 = zeros(6,6); % (6x6) position & velocity of LISA 1 N1 = zeros(6,1); % (6x1) M2 = zeros(6,6); %LISA 2 N2 = zeros(6,1); M3 = zeros(6,6); %LISA 3 N3 = zeros(6,1); case 3, %DSN & RELATIVE M = zeros(18,18); %LISA 1, 2 & 3 Postion & Velocity (18x18) N = zeros(18,1); %(18x1) otherwise, error('Invalid case (0 <= obs <= 3)'); end RMS1 = zeros(1); RMS2 = zeros(1); RMS3 = zeros(1); var1 = zeros(1); var2 = zeros(1); var3 = zeros(1); %---B. INNER LOOP for i = 1:length(ts), %loop over number of measurements for length of timespan %---READ IN TIME (1x1) ts(i); %each measurement corresponds to a time in span %---INTEGRATE STATE TRANSITION MATRIX, REFERENCE TRAJECTORY, & TRUTH 179 %TRAJECTORY FROM to->ti %NOTE: STM & REF. TRAJECTORY INTEGRATED ON THE NOMINAL ORBIT X* if i==1, Xstm = Xstmi'; %(1x600) no integration required for zero time change Sf = Xstm(end,:); %(1x600) Xstar=Xstaro; %(1x24) X = Xo; else [t, Xstm] = ode45('LISA_stmderive', [to ts(i)], Xstmi, odeset('RelTol', 1e-8,'AbsTol',1e-8)); [t, Xstar] = ode45('LISA_EOMderivs', [to ts(i)], Xstaro, odeset('RelTol', 1e-8,'AbsTol',1e-8)); [t, X] = ode45('LISA_EOMderivs', [to ts(i)], Xo, odeset('RelTol', 1e- 8,'AbsTol',1e-8)); Sf = Xstm(end,:); %(1x600) Xstar = Xstar(end, :); %1x24 Also, Xstar = Sf(1:24) X = X(end, :); %1x24 end %--Xstmi=Xstar plus the reshaped stm matrix as the I.C. for stmderive stm = reshape(Sf(25:600),24,24); %(24 x 24)the reshaped State Transition Matrix is columns 25:600 from Sf %---CALL RealObs.m & InterSCD.m to READ IN Measurement Deviations [Gi(i,:), L1tot(:,i), L2tot(:,i), L3tot(:,i)] = RealObs(to, ts(i), I, X, Xstar, S(i,1:3), S(i,4:6), S(i,7:9), S(i,10:12), S(i,19:21), S(i,22:24), B(1:3), 0); [GiR(i,:), L13tot(:,i), L21tot(:,i), L32tot(:,i)] = InterSCD(to, ts(i), I, X, Xstar, S(i,13:15), S(i,16:18), S(i,25:27), B(4:6), 0); %Total Measurement deviations switch(obs) case 0 %DSN & VLBI yi_tot(:,i) = [L1tot(:,i)]; %(4xi) yl_tot(:,i) = [L2tot(:,i)]; %(4xi) yq_tot(:,i) = [L3tot(:,i)]; %(4xi) case 1 %DSN, VLBI & RELATIVE yi_tot(:,i) = [L1tot(:,i); L13tot(:,i); L2tot(:,i); L21tot(:,i); L3tot(:,i); L32tot(:,i)]; %(18xi) case 2 %DSN only (range and range rate) yi_tot(:,i) = [L1tot(1); L1tot(2)]; %(2xi) yl_tot(:,i) = [L2tot(1); L2tot(2)]; %(2xi) yq_tot(:,i) = [L3tot(1); L3tot(2)]; %(2xi) case 3, %DSN & RELATIVE yi_tot(:,i) = [L1tot(1:2,i); L13tot(:,i); L2tot(1:2,i); L21tot(:,i); L3tot(1:2,i); L32tot(:,i)]; %(12xi) otherwise error('Invalid case (0 <= obs <= 3)'); end %---COMPUTE Hi_tilda & Hi %---H_tilda matrix, Observation Sensitivity (partial derivatives of the observations %over the State),or Observation-State Mapping Matrix, G = observations switch(obs) case 0, %DSN & VLBI [Hi_tilda(:,:,i), Hl_tilda(:,:,i), Hq_tilda(:,:,i)] = H_partials7(Gi(i,:), GiR(i,:), Xstar, 0); %CALCULATE THE OBSERVATION SENSITIVITY MATRIX OF EACH LISA Hi(:,:,i) = Hi_tilda(:,1:6,i)*stm(1:6,1:6); %(4x6) 180 Hl(:,:,i) = Hl_tilda(:,7:12,i)*stm(7:12,7:12); Hq(:,:,i) = Hq_tilda(:,13:18,i)*stm(13:18,13:18); %---READ IN Observation Noise Covariances (Ri) %construct as block diagonal matrices %---Weighting Matrix %NOTE: W = inv(Ri), Construct W matrices W1 = inv(diag(R1)); % (4x4) W2 = inv(diag(R2)); W3 = inv(diag(R3)); %return %---Scaled Weights W1sc = inv(diag(R1)); W2sc = inv(diag(R2)); W3sc = inv(diag(R3)); case 1, %DSN, VLBI, RELATIVE [Hi_tilda(:,:,i), Hl_tilda(:,:,i), Hq_tilda(:,:,i)] = H_partials7(Gi(i,:), GiR(i,:), Xstar, 1); Hi(:,:,i) = Hi_tilda(:,:,i)*stm(1:18,1:18); %(18x18) %---Weighting Matrix W1 = inv(diag(R1)); % (18x18) %---Scaled Weights W1sc = inv(diag(R1)); case 2, %DSN only (range and range rate) [Hi_tilda(:,:,i), Hl_tilda(:,:,i), Hq_tilda(:,:,i)] = H_partials7(Gi(i,:), GiR(i,:), Xstar, 2); %CALCULATE THE OBSERVATION SENSITIVITY MATRIX OF EACH LISA Hi(:,:,i) = Hi_tilda(:,1:6,i)*stm(1:6,1:6); %(2x6) Hl(:,:,i) = Hl_tilda(:,7:12,i)*stm(7:12,7:12); Hq(:,:,i) = Hq_tilda(:,13:18,i)*stm(13:18,13:18); %---Weighting Matrix W1 = inv(diag(R1)); % (4x4) W2 = inv(diag(R2)); W3 = inv(diag(R3)); %---Scaled Weights W1sc = inv(diag(R1)); W2sc = inv(diag(R2)); W3sc = inv(diag(R3)); case 3, %DSN & RELATIVE [Hi_tilda(:,:,i), Hl_tilda(:,:,i), Hq_tilda(:,:,i)] = H_partials7(Gi(i,:), GiR(i,:), Xstar, 3); Hi(:,:,i) = Hi_tilda(:,:,i)*stm(1:18,1:18); %(18x18) %---Weighting Matrix W1 = inv(diag(R1)); % (12x12) %---Scaled Weights W1sc = inv(diag(R1)); otherwise, error('Invalid case (0 <= obs <= 3)'); end %ACCUMULATE CURRENT OBSERVATION switch(obs) case 0, %DSN & VLBI SC=diag([1 1 1 1e-6 1e-6 1e-6]); %Scaling Matrix M1 = M1 + ((Hi(:,:,i)*SC)'*W1sc*(Hi(:,:,i)*SC)); %6x6, 2-D N1 = N1 + ((Hi(:,:,i)*SC)'*W1sc*yi_tot(:,i)); %6x1, 2-D 181 M2 = M2 + ((Hl(:,:,i)*SC)'*W2sc*(Hl(:,:,i)*SC)); N2 = N2 + ((Hl(:,:,i)*SC)'*W2sc*yl_tot(:,i)); M3 = M3 + ((Hq(:,:,i)*SC)'*W3sc*(Hq(:,:,i)*SC)); N3 = N3 + ((Hq(:,:,i)*SC)'*W3sc*yq_tot(:,i)); %*********************************************************** ** RMS1 = RMS1 + (yi_tot(:,i)'*W1sc*yi_tot(:,i)); %1x1, Weighted RMS RMS2 = RMS2 + (yl_tot(:,i)'*W2sc*yl_tot(:,i)); RMS3 = RMS3 + (yq_tot(:,i)'*W3sc*yq_tot(:,i)); RMSm = [RMS1 RMS2 RMS3]; var1 = var1 + yi_tot(:,i)'*W1sc*yi_tot(:,i); var2 = var2 + yl_tot(:,i)'*W2sc*yl_tot(:,i); var3 = var3 + yq_tot(:,i)'*W3sc*yq_tot(:,i); varm = [var1 var2 var3]; case 1 %DSN, VLBI & RELATIVE SC = diag([1 1 1 1e-6 1e-6 1e-6 1 1 1 1e-6 1e-6 1e-6... 1 1 1 1e-6 1e-6 1e-6]); M = M + ((Hi(:,:,i)*SC)'*(W1sc)*(Hi(:,:,i)*SC)); %18x18, 2-D N = N + ((Hi(:,:,i)*SC)'*(W1sc)*yi_tot(:,i)); %18x1, 2-D %*********************************************************** ** RMS1 = RMS1 + (yi_tot(1:6,i)'*(W1sc(1:6,1:6))*yi_tot(1:6,i)); %1x1, Weighted RMS RMS2 = RMS2 + (yi_tot(7:12,i)'*(W1sc(7:12,7:12))*yi_tot(7:12,i)); %1x1, Weighted RMS RMS3 = RMS3 + (yi_tot(13:18,i)'*(W1sc(13:18,13:18))*yi_tot(13:18,i)); %1x1, Weighted RMS RMSm = [RMS1 RMS2 RMS3]; var1 = var1 + yi_tot(1:6,i)'*(W1sc(1:6,1:6))*yi_tot(1:6,i); var2 = var2 + yi_tot(7:12,i)'*(W1sc(7:12,7:12))*yi_tot(7:12,i); var3 = var3 + yi_tot(13:18,i)'*(W1sc(13:18,13:18))*yi_tot(13:18,i); varm = [var1 var2 var3]; case 2 %DSN SC=diag([1 1 1 1e-6 1e-6 1e-6]); M1 = M1 + ((Hi(:,:,i)*SC)'*(W1sc)*(Hi(:,:,i)*SC)); %6x6, 2-D N1 = N1 + ((Hi(:,:,i)*SC)'*(W1sc)*yi_tot(:,i)); %6x1, 2-D M2 = M2 + ((Hl(:,:,i)*SC)'*(W2sc)*(Hl(:,:,i)*SC)); N2 = N2 + ((Hl(:,:,i)*SC)'*(W2sc)*yl_tot(:,i)); M3 = M3 + ((Hq(:,:,i)*SC)'*(W3sc)*(Hq(:,:,i)*SC)); N3 = N3 + ((Hq(:,:,i)*SC)'*(W3sc)*yq_tot(:,i)); %*********************************************************** ** RMS1 = RMS1 + (yi_tot(:,i)'*(W1sc)*yi_tot(:,i)); %1x1, Weighted RMS RMS2 = RMS2 + (yl_tot(:,i)'*(W1sc)*yl_tot(:,i)); RMS3 = RMS3 + (yq_tot(:,i)'*(W1sc)*yq_tot(:,i)); RMSm = [RMS1 RMS2 RMS3]; var1 = var1 + yi_tot(:,i)'*W1sc*yi_tot(:,i); var2 = var2 + yl_tot(:,i)'*W2sc*yl_tot(:,i); var3 = var3 + yq_tot(:,i)'*W3sc*yq_tot(:,i); varm = [var1 var2 var3]; case 3, %DSN & RELATIVE SC=diag([1 1 1 1e-6 1e-6 1e-6 1 1 1 1e-6 1e-6 1e-6... 1 1 1 1e-6 1e-6 1e-6]); M = M + ((Hi(:,:,i)*SC)'*(W1sc)*(Hi(:,:,i)*SC)); %18x18, 2-D N = N + ((Hi(:,:,i)*SC)'*(W1sc)*yi_tot(:,i)); %18x1, 2-D 182 %********************************************************** RMS1 = RMS1 + (yi_tot(1:4,i)'*(W1sc(1:4,1:4))*yi_tot(1:4,i)); %1x1, %Weighted RMS RMS2 = RMS2 + (yi_tot(5:8,i)'*(W1sc(5:8,5:8))*yi_tot(5:8,i)); %1x1, %Weighted RMS RMS3 = RMS3 + (yi_tot(9:12,i)'*(W1sc(9:12,9:12))*yi_tot(9:12,i)); %1x1, Weighted RMS RMSm = [RMS1 RMS2 RMS3]; var1 = var1 + yi_tot(1:4,i)'*W1sc(1:4,1:4)*yi_tot(1:4,i); var2 = var2 + yi_tot(5:8,i)'*W1sc(5:8,5:8)*yi_tot(5:8,i); var3 = var3 + yi_tot(9:12,i)'*W1sc(9:12,9:12)*yi_tot(9:12,i); varm = [var1 var2 var3]; otherwise error('Invalid case (0 <= obs <= 3)'); end if ts(i) < tf, if ts(i)== ts(1) to = ts(i); %time ti becomes to i = i + 1; Xstaro = Xstaro; %(1x24) the I.C. @ i=1 becomes the I.C. for next %reading Xstmi = Xstmi; %(600x1), STM becomes the I.C. for next reading Xo = Xo; %the I.C. @ i=1 becomes the I.C. for the next reading else to = ts(i); %time ti becomes to i = i + 1; Xstaro = Xstar; %(1x24) integrated trajectory becomes the I.C. for next %reading Xstmi = Sf'; %Sf' = (600x1), integrated STM becomes the I.C. for next %reading Xo = X; %integrated truth becomes the I.C. for the next reading end %---GO BACK TO THE BEGINNING OF B TO READ NEXT %OBSERVATION---% elseif ts(i) >= tf disp('---All measurements in batch have been read---'); end %PROCEED TO STEP C end %END OF INNER LOOP B, when all measurements in the batch are read disp('done with inner loop B'); %*********************************************************** ******* %OBSERVABILITY & MATRIX RANK TEST %IF FULL RANK & OBSERVABILITY IS GREATER THAN ZERO (i.e. not %singular %(det(M) is not equal to zero), to be singular, det(M)=0) then M is %invertible and observable %A test for the existence of exact solutions is rank(M) = rank([M N]), if %not equal then the answer is a LS solution, and is not exact %*********************************************************** ******* %---C. SOLVE THE NORMAL EQUATIONS %Alternative way to do the inversion of the Normal equations..... %Here is the straight up way, must pass the observability and matrix %rank test above to be able to invert %**************************** l = (length(ts)); %(1:length(ts)) %length of batch, i p = length(yi_tot(:,end)); % length of p-dimensional set of observation (usu. p> n switch(obs) 183 case 0 %DSN & VLBI rcond(M1); rcond(M2); rcond(M3); rank(M1); rank([M1 N1]); rank(M2); rank([M2 N2]); rank(M3); rank([M3 N3]); P1 = inv(M1); xhat1 = SC*P1*N1; P2 = inv(M2); xhat2 = SC*P2*N2; P3 = inv(M3); xhat3 = SC*P3*N3; n_xhat1 = norm(xhat1); n_xhat2 = norm(xhat2); n_xhat3 = norm(xhat3); n_xcorr = [n_xhat1 n_xhat2 n_xhat3]; xhat = [xhat1 xhat2 xhat3]; n_xhat = norm(xhat); %***RMS & Sample Variance RMSm; RMS = ([RMSm]./m).^(1/2); %(1x3) var = (1/(m-6)).*([varm]); %(1x3) case 1 %DSN, VLBI & RELATIVE rcond(M); rank(M); rank([M N]); P = inv(M); xhat = SC*P*N; n_xhat = norm(xhat); n_xhat1 = norm(xhat(1:6)); n_xhat2 = norm(xhat(7:12)); n_xhat3 = norm(xhat(13:18)); n_xcorr = [n_xhat1 n_xhat2 n_xhat3]; %***RMS & Sample Variance RMS = ([RMSm]./m).^(1/2); %(1x3) var = (1/(m-6)).*([varm]); %(1x3) var = (1/(m-18)).*([varm]); %(1x3) case 2 %DSN rcond(M1); rcond(M2); rcond(M3); rank(M1); rank([M1 N1]); rank(M2); rank([M2 N2]); rank(M3); rank([M3 N3]); P1 = inv(M1); xhat1 = SC*P1*N1; P2 = inv(M2); xhat2 = SC*P2*N2; P3 = inv(M3); xhat3 = SC*P3*N3; 184 n_xhat1 = norm(xhat1); n_xhat2 = norm(xhat2); n_xhat3 = norm(xhat3); n_xcorr = [n_xhat1 n_xhat2 n_xhat3]; xhat = [xhat1 xhat2 xhat3]; n_xhat = norm(xhat); %***RMS & Sample Variance RMS = ([RMSm]./m).^(1/2); %(1x3) var = (1/(m-6)).*([varm]); %(1x3) case 3, %DSN & RELATIVE rcond(M); rank(M); rank([M N]); P = inv(M); xhat = SC*P*N; n_xhat = norm(xhat); n_xhat1 = norm(xhat(1:6)); n_xhat2 = norm(xhat(7:12)); n_xhat3 = norm(xhat(13:18)); n_xcorr = [n_xhat1 n_xhat2 n_xhat3]; %***RMS & Sample Variance RMS = ([RMSm]./m).^(1/2); %(1x3) var = (1/(m-6)).*([varm]); %(1x3) var = (1/(m-18)).*([varm]); %(1x3) otherwise error('Invalid case (0 <= obs <= 3)') end %HAS PROCESS CONVERGED? niter = niter + 1; %dRMS = mean((abs(RMSold - RMS))); dRMS = mean((abs(RMSold - RMS))./RMSold); if dRMS >= tol, %THEN UPDATE REFERENCE STATE & ITERATE switch(obs) case 0 %DSN & VLBI Xhat1 = Xrefstart(1:6) + xhat1'; Xhat2 = Xrefstart(7:12) + xhat2'; Xhat3 = Xrefstart(13:18) + xhat3'; Xrefstart = [Xhat1 Xhat2 Xhat3 Xtru(1,19:24)]; %Updated reference/nominal trajectory for next iteration RMSold = RMS; case 1 %DSN, VLBI & RELATIVE Xhat = Xrefstart(1:18) + xhat'; Xrefstart = [Xhat Xtru(1,19:24)]; %Updated reference/nominal trajectory for next iteration RMSold = RMS; case 2 %DSN Xhat1 = Xrefstart(1:6) + xhat1'; Xhat2 = Xrefstart(7:12) + xhat2'; Xhat3 = Xrefstart(13:18) + xhat3'; Xrefstart = [Xhat1 Xhat2 Xhat3 Xtru(1,19:24)]; %Updated reference/nominal trajectory for next iteration RMSold = RMS; case 3, %DSN & RELATIVE Xhat = Xrefstart(1:18) + xhat'; Xrefstart = [Xhat Xtru(1,19:24)]; %Updated reference/nominal trajectory for next iteration RMSold = RMS; 185 otherwise error('Invalid case (0 <= obs <= 3)'); end %---GO BACK TO A = BEGINNING OF OUTER LOOP with new values Xstar and %xo_bar--% else dRMS < tol, %THE PROCESS HAS CONVERGED switch(obs) case 0 %DSN & VLBI Xhat1 = Xrefstart(1:6) + xhat1'; %updated reference state for LISA 1 Xhat2 = Xrefstart(7:12) + xhat2'; %for LISA 2 Xhat3 = Xrefstart(13:18) + xhat3'; %for LISA 3 Xup = [Xhat1 Xhat2 Xhat3]; %Updated reference/nominal trajectory Xest = [Xup Xtru(1,19:24)]; % give last Xhat as the resulting estimate from batch LS case 1 %DSN, VLBI & RELATIVE Xup = Xrefstart(1:18) + xhat'; %updated reference state for LISA 1,2,& 3 Xest = [Xup Xtru(1,19:24)]; %give last Xhat as the resulting estimate from batch LS case 2 %DSN Xhat1 = Xrefstart(1:6) + xhat1'; %updated reference state for LISA 1 Xhat2 = Xrefstart(7:12) + xhat2'; %for LISA 2 Xhat3 = Xrefstart(13:18) + xhat3'; %for LISA 3 Xup = [Xhat1 Xhat2 Xhat3]; %Updated reference/nominal trajectory Xest = [Xup Xtru(1,19:24)]; % give last Xhat as the resulting estimate from batch LS case 3, %DSN & RELATIVE Xup = Xrefstart(1:18) + xhat'; %updated reference state for LISA 1,2,& 3 Xest = [Xup Xtru(1,19:24)]; %give last Xhat as the resulting estimate from batch LS otherwise error('Invalid case (0 <= obs <= 3)'); end end %CAN END LOOPS %disp('done with outer loop A') end %END %---D. ERROR ANALYSIS switch(obs) case 0 %DSN & VLBI for j = 1:length(ts), e1(:,j) = yi_tot(:,j) - ([Hi(:,:,j)]*xhat1); %nxj best estimate of observation error %nx1, individually e2(:,j) = yl_tot(:,j) - ([Hl(:,:,j)]*xhat2); e3(:,j) = yq_tot(:,j) - ([Hq(:,:,j)]*xhat3); end obserr = [mean(e1,2) mean(e2,2) mean(e3,2)]; %nx3 %Unscaled Covariance SC=diag([1 1 1 1e-6 1e-6 1e-6]); %Scaling Matrix SC2 = (SC'*SC); P1_us = (SC2*P1); %(6x6) P2_us = (SC2*P2); P3_us = (SC2*P3); 186 Pcov = [diag(P1_us)' diag(P2_us)' diag(P3_us)']; %(1x18) sPcov = [sqrt(Pcov)]'; %Parameter Uncertainties PU_L1 =[sqrt(var(1)*abs(P1_us(1,1))); sqrt(var(1)*abs(P1_us(2,2))); sqrt(var(1)*abs(P1_us(3,3)));... sqrt(var(1)*abs(P1_us(4,4))); sqrt(var(1)*abs(P1_us(5,5))); sqrt(var(1)*abs(P1_us(6,6)))]; %(6x1) for LISA 1 Position & Velocity PU_L2 = [sqrt(var(2)*abs(P2_us(1,1))); sqrt(var(2)*abs(P2_us(2,2))); sqrt(var(2)*abs(P2_us(3,3)));... sqrt(var(2)*abs(P2_us(4,4))); sqrt(var(2)*abs(P2_us(5,5))); sqrt(var(2)*abs(P2_us(6,6)))]; %(6x1) for LISA 2 Position & Velocity PU_L3 = [sqrt(var(3)*abs(P3_us(1,1))); sqrt(var(3)*abs(P3_us(2,2))); sqrt(var(3)*abs(P3_us(3,3)));... sqrt(var(3)*abs(P3_us(4,4))); sqrt(var(3)*abs(P3_us(5,5))); sqrt(var(3)*abs(P3_us(6,6)))]; %(6x1) for LISA 3 Position & Velocity PU = [PU_L1 PU_L2 PU_L3]; %6x3 %***Position Errors err1r = norm(xhat1(1:3)); err2r = norm(xhat2(1:3)); err3r = norm(xhat3(1:3)); err = [err1r err2r err3r]; %***Velocity Errors err1v = norm(xhat1(4:6)); err2v = norm(xhat2(4:6)); err3v = norm(xhat3(4:6)); erv = [err1v err2v err3v]; case 1 %DSN, VLBI & RELATIVE for j = 1:length(ts), e1(:,j) = yi_tot(:,j) - ([Hi(:,:,j)]*xhat); %nxj best estimate of observation error %nx1, individually end obserr = [e1]; %Unscaled Covariances SC = diag([1 1 1 1e-6 1e-6 1e-6 1 1 1 1e-6 1e-6 1e-6... 1 1 1 1e-6 1e-6 1e-6]); SC2 = (SC'*SC); P_us = (SC2*P); %(18x18) Pcov = [diag(P_us)']; %Parameter Uncertainties PU_L1 =[sqrt(var(1)*abs(P_us(1,1))); sqrt(var(1)*abs(P_us(2,2))); sqrt(var(1)*abs(P_us(3,3)));... sqrt(var(1)*abs(P_us(4,4))); sqrt(var(1)*abs(P_us(5,5))); sqrt(var(1)*abs(P_us(6,6)))]; %(6x1) for LISA 1 Position & Velocity PU_L2 = [sqrt(var(2)*abs(P_us(7,7))); sqrt(var(2)*abs(P_us(8,8))); sqrt(var(2)*abs(P_us(9,9)));... sqrt(var(2)*abs(P_us(10,10))); sqrt(var(2)*abs(P_us(11,11))); sqrt(var(2)*abs(P_us(12,12)))]; %(6x1) for LISA 2 Position & Velocity PU_L3 = [sqrt(var(3)*abs(P_us(13,13))); sqrt(var(3)*abs(P_us(14,14))); sqrt(var(3)*abs(P_us(15,15)));... sqrt(var(3)*abs(P_us(16,16))); sqrt(var(3)*abs(P_us(17,17))); sqrt(var(3)*abs(P_us(18,18)))]; %(6x1) for LISA 3 Position & Velocity %***Position Errors err = [norm(xhat(1:3)) norm(xhat(7:9)) norm(xhat(13:15))]; %***Velocity Errors erv = [norm(xhat(4:6)) norm(xhat(10:12)) norm(xhat(16:18))]; case 2 %DSN for j = 1:length(ts), 187 e1(:,j) = yi_tot(:,j) - ([Hi(:,:,j)]*xhat1); %nxj best estimate of observation error %nx1, individually e2(:,j) = yl_tot(:,j) - ([Hl(:,:,j)]*xhat2); e3(:,j) = yq_tot(:,j) - ([Hq(:,:,j)]*xhat3); end obserr = [e1 e2 e3]; %Unscaled Covariance SC = diag([1 1 1 1e-6 1e-6 1e-6]); SC2 = (SC'*SC); P1_us = (SC2*P1); %(6x6) P2_us = (SC2*P2); P3_us = (SC2*P3); Pcov = [diag(P1_us)' diag(P2_us)' diag(P3_us)']; %Parameter Uncertainties PU_L1 =[sqrt(var(1)*abs(P1_us(1,1))); sqrt(var(1)*abs(P1_us(2,2))); sqrt(var(1)*abs(P1_us(3,3)));... sqrt(var(1)*abs(P1_us(4,4))); sqrt(var(1)*abs(P1_us(5,5))); sqrt(var(1)*abs(P1_us(6,6)))]; %(6x1) for LISA 1 Position & Velocity PU_L2 = [sqrt(var(2)*abs(P2_us(1,1))); sqrt(var(2)*abs(P2_us(2,2))); sqrt(var(2)*abs(P2_us(3,3)));... sqrt(var(2)*abs(P2_us(4,4))); sqrt(var(2)*abs(P2_us(5,5))); sqrt(var(2)*abs(P2_us(6,6)))]; %(6x1) for LISA 2 Position & Velocity PU_L3 = [sqrt(var(3)*abs(P3_us(1,1))); sqrt(var(3)*abs(P3_us(2,2))); sqrt(var(3)*abs(P3_us(3,3)));... sqrt(var(3)*abs(P3_us(4,4))); sqrt(var(3)*abs(P3_us(5,5))); sqrt(var(3)*abs(P3_us(6,6)))]; %(6x1) for LISA 3 Position & Velocity %***Position Errors err1r = norm(xhat1(1:3)); err2r = norm(xhat2(1:3)); err3r = norm(xhat3(1:3)); err = [err1r err2r err3r]; %***Velocity Errors err1v = norm(xhat1(4:6)); err2v = norm(xhat2(4:6)); err3v = norm(xhat3(4:6)); erv = [err1v err2v err3v]; case 3, %DSN & RELATIVE for j = 1:length(ts), e1(:,j) = yi_tot(:,j) - ([Hi(:,:,j)]*xhat); %nxj best estimate of observation error %nx1, individually end obserr = [e1]; %Unscaled Covariances SC = diag([1 1 1 1e-6 1e-6 1e-6 1 1 1 1e-6 1e-6 1e-6... 1 1 1 1e-6 1e-6 1e-6]); SC2 = (SC'*SC); P_us = (SC2*P); %(18x18) Pcov = [diag(P_us)']; %Parameter Uncertainties PU_L1 =[sqrt(var(1)*abs(P_us(1,1))); sqrt(var(1)*abs(P_us(2,2))); sqrt(var(1)*abs(P_us(3,3)));... sqrt(var(1)*abs(P_us(4,4))); sqrt(var(1)*abs(P_us(5,5))); sqrt(var(1)*abs(P_us(6,6)))]; %(6x1) for LISA 1 Position & Velocity PU_L2 = [sqrt(var(2)*abs(P_us(7,7))); sqrt(var(2)*abs(P_us(8,8))); sqrt(var(2)*abs(P_us(9,9)));... sqrt(var(2)*abs(P_us(10,10))); sqrt(var(2)*abs(P_us(11,11))); sqrt(var(2)*abs(P_us(12,12)))]; %(6x1) for LISA 2 Position & Velocity 188 PU_L3 = [sqrt(var(3)*abs(P_us(13,13))); sqrt(var(3)*abs(P_us(14,14))); sqrt(var(3)*abs(P_us(15,15)));... sqrt(var(3)*abs(P_us(16,16))); sqrt(var(3)*abs(P_us(17,17))); sqrt(var(3)*abs(P_us(18,18)))]; %(6x1) for LISA 3 Position & Velocity %***Position Errors err = [norm(xhat(1:3)) norm(xhat(7:9)) norm(xhat(13:15))]; %***Velocity Errors erv = [norm(xhat(4:6)) norm(xhat(10:12)) norm(xhat(16:18))]; otherwise error('Invalid case (0 <= obs <= 6)'); end %***ESTIMATE ERROR****************************************** fuzz = Xk - Xest; nfuzz = norm(fuzz); Xerror = Xest-Xtru; nXerror = norm(Xerror); nXest = norm(Xest); nXtru = norm(Xtru); %***ABSOLUTE OD ACCURACY************************************* DXr = [norm(Xest(1:3)-Xtru(1:3)), norm(Xest(7:9)-Xtru(7:9)), norm(Xest(13:15)-Xtru(13:15))]; %position DXv = [norm(Xest(4:6)-Xtru(4:6)), norm(Xest(10:12)-Xtru(10:12)), norm(Xest(16:18)-Xtru(16:18))]; %velocity ABS = [norm(DXr) norm(DXv)]; %***RELATIVE OD ACCURACY************************************** DL13=norm(Xest(1:3)-Xest(13:15))-norm(Xtru(1:3)-Xtru(13:15)); DL21=norm(Xest(1:3)-Xest(7:9))-norm(Xtru(1:3)-Xtru(7:9)); DL32=norm(Xest(7:9)-Xest(13:15))-norm(Xtru(7:9)-Xtru(13:15)); DL13v=norm(Xest(4:6)-Xest(16:18))-norm(Xtru(4:6)-Xtru(16:18)); DL21v=norm(Xest(10:12)-Xest(4:6))-norm(Xtru(10:12)-Xtru(4:6)); DL32v=norm(Xest(16:18)-Xest(10:12))-norm(Xtru(16:18)-Xtru(10:12)); DL = norm([DL13 DL21 DL32]); DLv = norm([DL13v DL21v DL32v]); REL = [DL DLv]; MCRUNbp.m (from Section 5.4 Orbit Accuracy) function [Mniter, MXest, Merr, Stderr, Merv, Stderv, MABSr, StdABSr, MABSv, StdABSv, MRELr, StdRELr, MRELv, StdRELv, MRMS, Mvar, MPU_L1, MPU_L2, MPU_L3, M_Pcov] = MCRUNbp() format long e x = 1; days = [1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90]; for d = [1 3 5 7 10 15 20 25 30 35 37 40 45 50 55 60 65 70 75 80 85 90]; d; x; %***SET THE APPROPRIATE START TIME TO THE %CASE************************* to_0 = 0; %EPOCH, X_truth(1,1); %start of timespan %***SET THE APPROPRIATE END TIME TO THE %CASE*************************** tf = 60*60*24*d; %days %***SET THE APPROPRIATE INTERVAL TIME TO THE %CASE********************** I = 60*60*24; %sec, skips every 24 hours for 1 day intervals 1 meas/day %***TIMESPAN********************************************** ts = to_0:I:tf; %timespan for INTERVALS 189 %*********************************************************** %***INITIAL CONDITION VECTOR FROM NomOrb_half.mat %*********************************************************** %Xo = [X_traj(1,2:7) X_traj(1, 8:13) X_traj(1, 14:19) X_traj(1,20:25)]; %1x24 Xo_0 = [9.071855731392364e-006, 1.481339600426210e+011, 2.475763859214967e+009,... -3.007344842423256e+004, 1.841210487873900e-012, 3.077216312703865e-014,... 2.487592427287221e+009, 1.503039437121128e+011, - 1.292020540364669e+009,... -2.963517041188058e+004, 2.476377264069181e+002, 4.268667617608362e+002,... -2.487592427287196e+009, 1.503039437121128e+011, - 1.292020540364669e+009,... -2.963517041188058e+004, -2.476377264069132e+002, - 4.268667617608363e+002,... -5.116553726769669e+010, 1.405761582972606e+011, 0,... -2.798843961236834e+004, -1.018695892245916e+004, 0]; %*********************************************************** %***REFERENCE START VECTOR WITH ERROR aka "FUZZ" %ADDED ON %*********************************************************** errorstart = 0*[-1 1 1 1e-6 1e-6 1e-6 -1 1 1 1e-6 1e-6 1e-6 -1 1 1 1e-6 1e- 6 1e-6].*randn(1,18); %NO Guess Error Xrefstart_0 = Xo_0 + [errorstart zeros(1,6)]; %*********************************************************** %READ IN SIGMA ERROR, BIASES & WEIGHTS FOR EACH CASE %*********************************************************** c = 299792458; %m/s speed of light tau = 1000; %sec, Doppler count time Astd = 1e-15; %Allan Standard deviation for a Doppler Count time =1000s pi1 = norm([Xrefstart_0(1:3) - Xrefstart_0(19:21)]); %m (1x1) range scalar for LISA 1 pi2 = norm([Xrefstart_0(7:9) - Xrefstart_0(19:21)]); pi3 = norm([Xrefstart_0(13:15) - Xrefstart_0(19:21)]); pi13 = norm([Xrefstart_0(1:3) - Xrefstart_0(13:15)]); %m (1x1) range scalar for LISA Leg 13 pi21 = norm([Xrefstart_0(7:9) - Xrefstart_0(1:3)]); pi32 = norm([Xrefstart_0(13:15) - Xrefstart_0(7:9)]); RTLT1 = (2*pi1)/c; %round-trip light time for LISA 1 RTLT2 = (2*pi2)/c; RTLT3 = (2*pi3)/c; RTLT13 = (2*pi13)/c; %round-trip light time for LISA Leg 13 RTLT21 = (2*pi21)/c; RTLT32 = (2*pi32)/c; %range error due to clock instability ci1 = ((sqrt(2))*c*RTLT1*Astd)*sqrt(60/(1000)); ci2 = ((sqrt(2))*c*RTLT2*Astd)*sqrt(60/(1000)); ci3 = ((sqrt(2))*c*RTLT3*Astd)*sqrt(60/(1000)); ci = norm([ci1 ci2 ci3]); ci13 = ((sqrt(2))*c*RTLT13*Astd)*sqrt(60/(1000)); ci21 = ((sqrt(2))*c*RTLT21*Astd)*sqrt(60/(1000)); ci32 = ((sqrt(2))*c*RTLT32*Astd)*sqrt(60/(1000)); cileg = norm([ci13 ci21 ci32]); M1 = RTLT1/tau; %constant M M2 = RTLT2/tau; M3 = RTLT3/tau; M13 = RTLT13/tau; %constant M M21 = RTLT21/tau; M32 = RTLT32/tau; %range rate frequency instability to range rate error f1 = ((sqrt(2+ log2(M1)))*c*Astd)*sqrt(60/(1000)); f2 = ((sqrt(2+ log2(M2)))*c*Astd)*sqrt(60/(1000)); f3 = ((sqrt(2+ log2(M3)))*c*Astd)*sqrt(60/(1000)); 190 f = norm([f1 f2 f3]); %***SIGMA ERRORS, S sr = .60*sqrt(60/(30*60)); %m, for 60-sec average srr = 3e-5*sqrt(60/(30*60)); %m/s, range rate error for 60-sec average sa = 5e-8*sqrt((30*60)/(30*60)); %rad, right ascension error for 30-min average sd = 5e-8*sqrt((30*60)/(30*60)); %rad, declination error for 30-min average sir = .60*sqrt(60/(30*60));%3*sqrt(1/60); %m, interspacecraft range error for 1-sec average sirr = 3e-5*sqrt(60/(30*60));%1.5e-4*sqrt(1/60); %m/s, interspacecraft range rate error for 1-sec average %*********************************************************** %Observation Noise Covariances, Ri's %*********************************************************** %With Noise************************************* R0a = [(sr + ci1) (srr + f1) sd sa].^2; R0b = [(sr + ci2) (srr + f2) sd sa].^2; R0c = [(sr + ci3) (srr + f3) sd sa].^2; R1 = [(sr+ ci1) (srr + f1) sd sa (sir + ci13) (sirr) (sr+ ci2) (srr+ f2)... sd sa (sir+ci21) (sirr) (sr+ ci3) (srr+ f3) sd sa (sir+ci32) (sirr)].^2; R2a = [(sr+ ci1) (srr+ f1)].^2; R2b = [(sr+ ci2) (srr+ f2)].^2; R2c = [(sr+ ci3) (srr+ f3)].^2; R3 = [(sr+ ci1) (srr+ f1) (sir + ci13) (sirr) (sr+ ci2) (srr+ f2)... (sir + ci21) (sirr) (sr+ ci3) (srr+ f3) (sir + ci32) (sirr)].^2; %*********************************************************** %BIASES %*********************************************************** %RANGE RATE pdot1i = ((((Xrefstart_0(1) - Xrefstart_0(19))*(Xrefstart_0(4) - Xrefstart_0(22))) + ((Xrefstart_0(2) - Xrefstart_0(20))*... (Xrefstart_0(5) - Xrefstart_0(23))) + ((Xrefstart_0(3) - Xrefstart_0(21))*(Xrefstart_0(6) - Xrefstart_0(24))))/pi1); %m/s pdot2i =((((Xrefstart_0(7) - Xrefstart_0(19))*(Xrefstart_0(10) - Xrefstart_0(22))) + ((Xrefstart_0(8) - Xrefstart_0(20))*... (Xrefstart_0(11) - Xrefstart_0(23))) + ((Xrefstart_0(9) - Xrefstart_0(21))*(Xrefstart_0(12) - Xrefstart_0(24))))/pi2); pdot3i =((((Xrefstart_0(13) - Xrefstart_0(19))*(Xrefstart_0(16) - Xrefstart_0(22))) + ((Xrefstart_0(14) - Xrefstart_0(20))*... (Xrefstart_0(17) - Xrefstart_0(23))) + ((Xrefstart_0(15) - Xrefstart_0(21))*(Xrefstart_0(18) - Xrefstart_0(24))))/pi3); pdot13i = ((((Xrefstart_0(1) - Xrefstart_0(13))*(Xrefstart_0(4) - Xrefstart_0(16))) + ((Xrefstart_0(2) - Xrefstart_0(14))*... (Xrefstart_0(5) - Xrefstart_0(17))) + ((Xrefstart_0(3) - Xrefstart_0(15))*(Xrefstart_0(6) - Xrefstart_0(18))))/pi13); %m/s pdot21i =((((Xrefstart_0(7) - Xrefstart_0(1))*(Xrefstart_0(10) - Xrefstart_0(4))) + ((Xrefstart_0(8) - Xrefstart_0(2))*... (Xrefstart_0(11) - Xrefstart_0(5))) + ((Xrefstart_0(9) - Xrefstart_0(3))*(Xrefstart_0(12) - Xrefstart_0(6))))/pi21); pdot32i =((((Xrefstart_0(13) - Xrefstart_0(7))*(Xrefstart_0(16) - Xrefstart_0(10))) + ((Xrefstart_0(14) - Xrefstart_0(8))*... (Xrefstart_0(17) - Xrefstart_0(11))) + ((Xrefstart_0(15) - Xrefstart_0(9))*(Xrefstart_0(18) - Xrefstart_0(12))))/pi32); %range error due to unknown time-tag error dT = 1e-6; %Station clock epoch error g1 = dT*pdot1i; g2 = dT*pdot2i; g3 = dT*pdot3i; g13 = dT*pdot13i; g21 = dT*pdot21i; g32 = dT*pdot32i; 191 %clock rate offset rT = 5e-14; %Station clock rate offset h1 = c*RTLT1*rT; h2 = c*RTLT2*rT; h3 = c*RTLT3*rT; h13 = c*RTLT13*rT; h21 = c*RTLT21*rT; h32 = c*RTLT32*rT; br = 2.05; %2 m instrument bias + 5 cm total bias bir = 2; %10; %10 m bias for each LISA B = [(br+g1+h1) (br+g2+h2) (br+g3+h3) (bir + g13 + h13) (bir + g21 + h21) (bir + g32 + h32)]; %biases 1x6 %B = [0 0 0 0 0 0]; %No biases 1x6 %*********************************************************** %***READ IN OTHER INPUTS %*********************************************************** plt = 0; %*********************************************************** %***CALL LISA_BPtot4.m to RUN DIFFERENT CASES %*********************************************************** %for i = 1:35, %100 random runs of each case for i = 1, %1 run of each for full noise, 100 random runs of each case for k = 1:length(ts) %S = ([sr sr sr srr srr srr sd sd sd 0 0 0 sir sir sir sirr sirr sirr ci1 ci2 ci3 f1 f2 f3 ci13 ci21 ci32]'*zeros(1, k+1))'; %No Noise, nx27 S = ([sr sr sr srr srr srr sd sd sd sa sa sa sir sir sir sirr sirr sirr ci1 ci2 ci3 f1 f2 f3 ci13 ci21 ci32]'*ones(1, k+1))'; %Full Noise, nx27 %S = randn(k+1, 27).*([sr sr sr srr srr srr sd sd sd sa sa sa sir sir sir sirr sirr sirr ci1 ci2 ci3 f1 f2 f3 ci13 ci21 ci32]'*ones(1, k+1))'; %Random Noise, nx27 %errorstart = randn(k+1, 18).*([-1 1 1 1e-6 1e-6 1e-6 -1 1 1 1e-6 1e- 6 1e-6 -1 1 1 1e-6 1e-6 1e-6]'*(1000*ones(1,k+1))); %NO Guess Error %errorstart = randn(1, 18).*([-1 1 1 1e-6 1e-6 1e-6 -1 1 1 1e-6 1e-6 1e-6 -1 1 1 1e-6 1e-6 1e-6]*(1000.*ones(1,1))) Xrefstart_0 = Xo_0 + [errorstart zeros(1,6)]; end %CASE 0: DSN & VLBI (All absolute data) %[Xest, ABS, REL, RMS, var, PU_L1, PU_L2, PU_L3, err, erv, niter, Pcov] = MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R1, R2, R3, S, B, plt, obs) [Xest0(i,:), ABS0(i,:), REL0(i,:), RMS0(i,:), var0(i,:),PU0_L1(:,i), PU0_L2(:,i), PU0_L3(:,i), err0(i,:), erv0(i,:), niter0(i,:), Pcov0(i,:)] = MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R0a, R0b, R0c, S, B, plt, 0); %CASE 1: DSN, VLBI & RELATIVE [Xest1(i,:), ABS1(i,:), REL1(i,:), RMS1(i,:),var1(i,:),PU1_L1(:,i), PU1_L2(:,i), PU1_L3(:,i), err1(i,:), erv1(i,:), niter1(i,:), Pcov1(i,:)] = MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R1, R1, R1, S, B, plt, 1); %CASE 2: DSN [Xest2(i,:), ABS2(i,:), REL2(i,:), RMS2(i,:), var2(i,:), PU2_L1(:,i), PU2_L2(:,i), PU2_L3(:,i), err2(i,:), erv2(i,:), niter2(i,:), Pcov2(i,:)] = MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R2a, R2b, R2c, S, B, plt, 2); %CASE 3: DSN & RELATIVE [Xest3(i,:), ABS3(i,:), REL3(i,:), RMS3(i,:), var3(i,:), PU3_L1(:,i),PU3_L2(:,i), PU3_L3(:,i), err3(i,:), erv3(i,:), niter3(i,:), Pcov3(i,:)] = MCLISA_bp(to_0, tf, I, Xrefstart_0, Xo_0, R3, R3, R3, S, B, plt, 3); end %*********************************************************** %****ERROR ACCUMULATION & ANALYSIS %*********************************************************** for j = x, 192 %4 Case [case 0(DSN/VLBI),case 1(DSN/VLBI/REL),case 2(DSN),case 3(DSN/REL)] niter_tot = [niter0(:,:), niter1(:,:), niter2(:,:), niter3(:,:)]; %nx4 M_niter(j,:) = mean(niter_tot, 1); %1x6 Xest = [Xest0(:,:), Xest1(:,:), Xest2(:,:), Xest3(:,:)]; %n x 144 (1x24 each) M_Xest(j,:) = mean(Xest,1); %1x144 ABSr = [ABS0(:,1), ABS1(:,1), ABS2(:,1), ABS3(:,1)]; %ix6 M_ABSr(j,:) = mean(ABSr,1); %1x6 std_ABSr(j,:) = std(ABSr, 0, 1); %1x6 ABSv = [ABS0(:,2), ABS1(:,2), ABS2(:,2), ABS3(:,2)]; %ix6 M_ABSv(j,:) = mean(ABSv,1); %1x6 std_ABSv(j,:) = std(ABSv, 0, 1); %1x6 RELr = [REL0(:,1), REL1(:,1), REL2(:,1), REL3(:,1)]; %ix6 M_RELr(j,:) = mean(RELr, 1); %1x6 std_RELr(j,:) = std(RELr, 0, 1); %1x6 RELv = [REL0(:,2), REL1(:,2), REL2(:,2), REL3(:,2)]; %ix6 M_RELv(j,:) = mean(RELv,1); %1x6 std_RELv(j,:)= std(RELv, 0, 1); %1x6 RMS = [RMS0(:,:), RMS1(:,:), RMS2(:,:), RMS3(:,:)]; %ix18 (1x3 for each case) M_RMS(j,:) = mean(RMS,1); %1x18 var = [var0(:,:), var1(:,:), var2(:,:), var3(:,:)]; %ix18 (1x3 for each case) M_var(j,:) = mean(var,1); %1x18 PU_L1 = [PU0_L1(:,1);... PU1_L1(:,1);... PU2_L1(:,1);... PU3_L1(:,1)]; %ix144(6x1 each) M_PU_L1(:,j) = mean(PU_L1,2); PU_L2 = [PU0_L2(:,1);... PU1_L2(:,1);... PU2_L2(:,1);... PU3_L2(:,1)]; M_PU_L2(:,j) = mean(PU_L2,2); PU_L3 = [PU0_L3(:,1);... PU1_L3(:,1);... PU2_L3(:,1);... PU3_L3(:,1)]; M_PU_L3(:,j) = mean(PU_L3,2); err = [err0(:,:), err1(:,:), err2(:,:), err3(:,:)]; %ix18 (1x3 each) M_err(j,:) = mean(err,1); %1x18 std_err(j,:) = std(err, 0, 1); erv = [erv0(:,:), erv1(:,:), erv2(:,:), erv3(:,:)]; %ix18 M_erv(j,:) = mean(erv,1); %1x18 std_erv(j,:) = std(erv, 0, 1); Pcov = [Pcov0(1,:), Pcov1(1,:), Pcov2(1,:), Pcov3(1,:)]; %ix(18x3)= ix54 M_Pcov(j,:) = mean(Pcov,1); %1x(18x4) if x < length(days) x = x + 1 else x = length(days) end end Mniter = M_niter(:,:); MXest = M_Xest(:,:); Merr = M_err(:,:); Stderr = std_err(:,:); 193 Merv = M_erv(:,:); Stderv = std_erv(:,:); MABSr = M_ABSr(:,:); StdABSr = std_ABSr(:,:); MABSv = M_ABSv(:,:); StdABSv = std_ABSv(:,:); MRELr = M_RELr(:,:); StdRELr = std_RELr(:,:); MRELv = M_RELv(:,:); StdRELv = std_RELv(:,:); MRMS = M_RMS(:,:); Mvar = M_var(:,:); MPU_L1 = M_PU_L1(:,:); MPU_L2 = M_PU_L2(:,:); MPU_L3 = M_PU_L3(:,:); M_Pcov = M_Pcov(:,:); end 194 Bibliography Listed in order of appearance in text [1] LISA NASA/JPL Homepage. URL: http://lisa.jpl.nasa.gov [2] Sorenson, H. W., ?Least-squares estimation: from Gauss to Kalman,? IEEE Spectrum, Vol. 7, July 1970, pp. 63-68. [3] Folkner, W.M., "The Laser Interferometer Space Antenna Mission," AIAA Space 2000 Conference and Exposition, Long Beach, CA, 19-21 September 2000 (AIAA 2000-5129). [4] Hechler, F. and Folkner, W. M., ?Mission Analysis for the Laser Interferometer Space Antenna (LISA) Mission,? Advanced Space Research, Vol. 32, No.7, Elsevier Ltd., Great Britain, 2003, pp. 1277-1282. [5] Vallado, David A., and Alfano, Salvatore, ?A Future Look at Space Surveillance and Operations,? Space Surveillance Workshop, Washington, DC, October 20-23, 1998. [6] Ziemer, John K., Merkowitz, Stephen M., ?Microthrust Propulsion for the LISA Mission,? The 40 th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit , Fort Lauderdale, FL., 11-14 July 2004 (AIAA 2004-3439). [7] Sweetser, Theodore, H., ?An End-to-end Trajectory Description of the LISA Mission,? Classical and Quantum Gravity 22, Institute of Physics Publishing Ltd., United Kingdom, April 2005, pp. 429-435. [8] Sweetser, Theodore, H., "LISA Mission Description," Internal Memo NASA/JPL, Version 2.0., Oct. 28, 2005. [9] LISA newsletter Vol. 1, No. 1, LISA NASA/JPL Homepage, URL: http://www.lisa.nasa.gov/newsletter/newsletter_200404.pdf, April 2004 [10] Irion, R., ?Gravitational Wave Hunters Take Aim at the Sky.? Science, Vol. 297, August 16, 2002, pp. 1113-1115. 195 [11] LISA Science Brochure, LISA NASA/JPL Homepage, URL: http://www.lisa.nasa.gov/LISAbrochure.pdf [12] Webster, Guy, (11/20/2001). NASA Spacecraft to Hunt for Elusive Gravity Ripples, Cassini-Huygens News Releases, URL: http://saturn.jpl.nasa.gov/news/press- releases-01/20011120-pr-a.cfm [13] GRACE: Gravity Recovery and Climate Experiment (2/11/2004) GRACE Gravity Measurement, URL: http://www.csr.utexas.edu/grace/science/gravity_measurement.html [14] Hughes, Steven P., ?Preliminary Optimal Orbit Design for the Laser Interferometer Space Antenna (LISA),? 25 th Annual AAS Guidance and Control Conference, Breckenridge CO, Feb. 2002. [15] Bender. P.L., ?LISA: A proposed joint ESA-NASA gravitational-wave mission?: Gravitational Waves, Institute of Physics Publishing Ltd., Philadelphia, 2001, Chapter 10, pp. 115-151. [16] Laser Interferometer Space Antenna: A cornerstone Mission for the Observation of Gravitational Waves, Systems and Technology Study (STS) Report, 11 July, 2000, URL: ftp://ftp.rzg.mpg.de/pub/grav/lisa/sts/sts_1.05.pdf [17] Fleck, Melissa E. and Starin, Scott R., ?Evaluation of a Drag-Free Control Concept for Missions in Low Earth Orbit,? AIAA Guidance, Navigation and Control Conference and Exhibit, Austin, Texas, 11-14 August 2003, (AIAA 2003-5747). [18] Dhurandhar, S. V., Rajesh, N.K., Kosti, S, and Vinet, J-Y, "Fundamentals of the LISA Stable Flight Formation," Classical and Quantum Gravity, Vol. 22, Institute of Physics Publishing Ltd., United Kingdom, Jan. 2005, pp. 481- 487. [19] Folkner, W. M., et al., ?LISA orbit selection and stability.? Classical and Quantum Gravity, Vol. 14, Institute of Physics Publishing Ltd., United Kingdom, 1997, pp.1405-1410. [20] LISA Mission Report from 35 th COSPAR meeting, found at LISA ESA Homepage, URL: http://lisa.esa.int/science-e/www/area/index.cfm?fareaid=27, pp.113- 116. 196 [21] NASA/JPL Space Technology 7 (ST7) About ST7 Overview 5/2004 URL: http://nmp.jpl.nasa.gov/st7/ABOUT/About_index.html. [22] ESA: Science and Technology Webpage LISA-Pathfinder Homepage Factsheet: 5/31/2005, URL: http://sci.esa.int/science-e/www/object/index.cfm?fobjectid=31714 -Instruments: last update 10/21/2004, URL: http://sci.esa.int/science-e/www/object/index.cfm?fobjectid=35589 -Summary: 5/31/2005, URL: http://sci.esa.int/sciencee/www/object/index.cfm?fobjectid=31431 [23] NASA Beyond Einstein: Great Observatories Homepage URL: http://universe.nasa.gov/program/observatories.html [24] NASA: GSFC Constellation X Homepage last update 1/6/2006 URL: http://constellation.gsfc.nasa.gov/ (checked 3/25/2005) [25] The TRIAD Spacecraft last updated 11/25/2001 URL: http://www- spof.gsfc.nasa.gov/Education/wtriad.html [26] NASA/JPL SIM Planetquest Mission Homepage URL: http://planetquest.jpl.nasa.gov/SIM/sim_index.html (checked 3/25/2005) [27] Savage, Don, et. al. NASA Gravity Probe B mission enters science phase, ready to test Einstein's theory NASA News Release 9/7/2004 URL: http://www.nasa.gov/centers/marshall/news/news/releases/2004/04-228.html [28] NASA Gravity Probe B Homepage URL: http://www.gravityprobeb.com/ [29] GFZ Potsdam Department 1, CHAMP Mission Homepage 11/29/2005 URL: http://www.gfz-potsdam.de/pb1/op/champ/main_CHAMP.html (checked 3/25/2005) [30] NASA/JPL Cassini-Huygens Homepage Overview: URL: http://saturn.jpl.nasa.gov/overview/index.cfm [31] Danzmann, K., & Rudiger, A., "LISA technology- concept, status, prospects," 197 Classical and Quantum Gravity, Vol. 14, Institute of Physics Publishing Ltd., United Kingdom, 2003, pp. S1-S9. [32] Crowder, J. & Corning, N. J., "Beyond LISA: Exploring Future Gravitational Wave Missions," 17 October 2005, URL: http://www.arxiv.org/PS_cache/gr- qc/pdf/0506/0506015.pdf. [33] Wiesel, William, E., Spaceflight Dynamics 2 nd ed., Irwin/McGraw-Hill Companies Inc., New York, NY, 1997. [34] Tapley, Byron D. ?Statistical Orbit Determination Theory,? Recent Advances in Dynamical Astronomy, D. Reidel Publishing Company, Boston, USA, 1973, pp. 396- 425. [35] Vallado, David A., Fundamentals of Astrodynamics and Applications, Space Technology Series ed. Wiley J. Larson, The McGraw-Hill Companies, Inc., New York, NY, 1997. [36] Battin, Richard, H., An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition, American Institute of Aeronautics and Astronautics, Inc., Reston, VA, 1999. [37] Chobotov, Vladimir, A. editor, Orbital Mechanics: AIAA Education Series, American Institute of Aeronautics and Astronautics, Inc., Washington, DC, 1991. [38] Boulet, Dan L. Jr., Methods of Orbit Determination for the Micro Computer, Willmann-Bell, Inc., Richmond, VA., 1991. [39] Vallado, David A., and Alfano, Salvatore, ?A Future Look at Space Surveillance and Operations,? Space Surveillance Workshop, Washington, DC, October 20-23, 1998. [40] Montenbruck, O., and Gill, E., Satellite Orbits- Models, Methods, and Applications, Springer-Verlag, Berlin, Germany, 2000. [41] Thornton, Catherine L., Border, James S., ?Radiometric Tracking Techniques for Deep Space Navigation?, Jet Propulsion Laboratory, California Institute of Technology, Oct. 2000. 198 [42] Haeberle, Dennis W., Spencer, David B., and Ely, Todd A.,?Preliminary Investigation of Interplanetary Navigation Using Large Antenna Arrays with Varied Baselines.? AIAA/AAS Astrodynamics Specialist Conference and Exhibit: Providence, RI., 16-19 August 2004, (AIAA 2004-4744) . [43] JPL:Radio Science Systems Group Website: ?Basics of Spaceflight,? http://www2.jpl.nasa.gov/basics/. -JPL Radio Astronomy DSN [JPL: Basics of Spaceflight ?Telecommunications? http://www2.jpl.nasa.gov/basics/bsf10-1.html [44] Portock, Brian, M., Baird, Darren, T., et. al., ?Mars Exploration Rovers Cruise Orbit Determination,? (AIAA 2004-4981), AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Providence, RI, 16-19 August 2004 (AIAA 2004-4981). [45] Tapley, Byron D., ?Fundamentals of Orbit Determination,? Lecture Notes in Earth Sciences, Vol. 25: Theory of Satellite Geodesy and Gravity Field Determination ed. F. Sanso and R. Rummel, Springer-Verlag, New York, NY, 1989, pp. 235-260. [46] Noton, Maxwell, Spacecraft Navigation and Guidance, Advances in Industrial Control Series, Springer-Verlag London Limited, Great Britain, 1998. [47] Tapley, Byron D., Schutz, Bob E., Born, George H., Statistical Orbit Determination, Elsevier Academic Press, New York, NY, 2004. [48] Bate, Roger, R., Mueller, Donald, D., and White, Jerry, E., Fundamentals of Astrodynamics, Dover Publications, Inc. New York, NY, 1971. [49] Gelb, A. ed. and Technical Staff, The Analytic Sciences Corporation, Applied Optimal Estimation, The M.I.T. Press, MIT Cambridge, MA, 1974. [50] Devore, Jay, L., Prabability and Statistics: For Engineering and the Sciences, 6 th ed., Brooks/Cole Thomson Learning, Inc., Belmont, CA, 2004. 199