Satellite Data Retrieval, Reference Frames, Numerical and Analytical Orbital Propagation SSD Individual Assignment RMIT University Figure 1: Example errors between a ground truth ephemeris and predicted (starlink) orbit. Introduction Often, path finding and feasibility projects for satellite systems will entail identifying a case study set of satellites, constellations, formations, etc. These case study satellites will be what the hypothetical control system, sensor, telecom hardware will evaluated for performance on. What makes a good case study satellite varies based upon what problem we wish to evaluate the performance of, but usually, matching orbital regime, period, nodal precession rate, and exposure to certain space environment factors is essential. Once an orbital case study has been identified, the next logical step is obtaining orbital data to actually use in simulation. Extreme precision and accuracy is not always essential, especially for very short time period studies. For example, an encounter with a hypothetical debris field may only last a few seconds each time it occurs. If we need only consider the first encounter, a simple central potential + J2 perturbation is likely satisfactory. If however, we want to identify the likely directions from which RSOs will approach our satellite over the course of a week, accurate propagators, significantly more orbital data for our case study satellite (as well as all other satellites in a similar regime) is essential. With this in mind, this assignment will evaluate (and build) your ability to conduct such a case study. You will use a common source of satellite orbital data (Space-Track.org) to obtain historical NORAD Two Line Element (TLE) data, and compare the resultant orbital state vectors (produced by SGP-4) against Ephemeris data provided by Space-X, as well as against increasingly more simple orbital models. The errors between these models will be compared in the body frame of a ground-truth orbit, and a fundamental understanding of orbital mechanics will be required in order to explain the distributions of error along the basis vectors of the body frame. Additionally, this assignment will require an understanding of the arithmetic behind these operations, going beyond the ability to use software. You will need to know fundamental equations that convert 1 between orbital state vectors and Keplerian elements, conduct J2 perturbation calculations, and calculate body frame transformation matrices. 1 Ephemeris Data Inside the ephemeris data you have been provided, the following columns are present: Table 1: Example Ephemeris Data yr yr_day r_x r_y r_z v_x v_y v_z 2023 149.4824 -5188.671 2719.162 -3702.620 0.005249263 -6.104685 -4.496573 The actual data you have has more digits of precision. Note that instead of MM/DD/YYYY HH:MM:SS, time steps are provided in day of year. This is fine for analysis, however should you wish to convert this to a more human readable format, a code snippet is provided in section 4; you may convert it to your language of choice (and realistically, your language of choice will probably have a built in function to do this). Read the 18th Space Control Squadronss book SPACEFLIGHT SAFETY HANDBOOK FOR SATELLITE OPERATORS overview on ephemeris formats. 1. Based on this, what is the coordinate frame in which the data you have been provided measured in? (1 pt) 2. Is it inertial, quasi-inertial, or earth fixed? (1 pt) Question 1 The ephemeris data you have been provided is not directly downloadable from Space-Track.org it is concatenated from several ephemeris data-sets downloaded over multiple days. There are some discontinuities between these concatenated data-sets, where the end of one ephemeris set is not perfectly lined up with the start of the next. This is particularly visible if you look at the magnitude of velocity calculated by taking differences of positions at different time steps, and further taking the magnitude of the differences of those velocities at different time steps (diff(...) and diff(diff(...)), respectively, in MATLAB or Python). Note that this is different to velocity provided in the ephemeris data. An example of this is provided below: (a) (b) Figure 2: Examples of discontinuities between concatenated ephemeris data-sets. 2 In your ephemeris data-set: 1. Determine the time-step size (1 pt) 2. Plot the numerically calculated ?r?? as described in the process above (4 pts) 3. Plot the numerically calculated ?r? as described in the process above (4 pts) 4. Are there discontinuities (look very carefully)? If so, highlight them on each plot, and explain what might be the cause. (2 pts) 5. Based on the average altitude the satellite is operating at, do a rough calculation (and show your working out!) of expected gravitational acceleration at that altitude and velocity of a satellite in a perfect circular orbit with semi-major axis equal to that altitude; does these match up with what you plotted? (4 pt) Question 2 2 TLE Data and SGP4 Propagation Error As mentioned before, TLE data for a given satellite can be found on Space-Track. In the ephemeris data-set you have been provided, the name of the file is the NORAD ID of the Starlink satellite youve been assigned. Create a Space-Track account, and read up on their api that is used for accessing TLE data. The API is essentially formulating a URL that tells the website what data you want to access. To match the epoch of the provided ephemeris data, well need to get a historical TLE. Most satellites have a TLE released two or three times a day, and it has been a while since the 150th day of year stated in the ephemeris data youve been provided. To get historical data, you need to use gp_history. Additionally, given that we want to use SGP4, well need two line elements (TLE), well need to filter by NORAD ID, and sort by TLE age. Additionally, to prevent Space-Track throwing an error, well need to set a limit on the number of total TLEs to download. This will build up our query URL. As an example, I am considering the starlink satellite with NORAD ID 45659 the first part of the query URL is then https://www.space-track.org/basicspacedata/query/class/gp_history/format/tle/ NORAD_CAT_ID/45659 Then, to sort by epoch (now being at the top) I can orderby with the argument EPOCH desc: /orderby/EPOCH desc and then set a limit of 3 (an arbitrary number) /limit/3 Once youve made a Space-Track account try concatenating these urls, and you should get something like this, using this query url: 1 45659U 20035C 23194.75645042 .00005069 00000-0 35883-3 0 9991 2 45659 53.0537 100.8186 0001405 83.0845 277.0304 15.06410472172133 1 45659U 20035C 23193.76125188 .00021154 00000-0 14335-2 0 9996 2 45659 53.0534 105.2863 0001258 88.8289 271.2844 15.06433571171984 1 45659U 20035C 23193.16413945 .00009987 00000-0 68827-3 0 9997 2 45659 53.0530 107.9661 0001514 99.0773 261.0387 15.06405451171897 Youll need to make sure you are signed into your account in another tab before clicking that link, otherwise the cookies Space-Track uses to know which account a request is coming from wont exist. The actual values you get wont be the same as those above, as NOW is probably mid-semester for you, but I am writing this part of the assignment before the semester has started! For the satellite 45659, the first epoch of the ephemeris data-set is 3 2023 149.48243055555557 To predict the location of the satellite at that point in time using SGP4, a TLE with a start epoch as close to (but not later than!) that epoch is required. For this satellite, that TLE is 1 45659U 20035C 23148.71170837 .00019754 00000-0 13406-2 0 9993 2 45659 53.0537 307.5228 0001108 92.3053 267.8062 15.06420255165590 Using the query URL formatting method above, extract three TLEs who epochs are closest to (but not later than!) the start epoch of the ephemeris data you have been provided. If there are repeated epochs, be sure to ignore them! (3 pts) Question 3 TLEs put into SGP4 work in the true equator, mean equinox coordinate frame. However, true equator has two definitions: True Equator of Date in this option, the epoch of the TEME frame is always the same as the epoch of the associated ephemeris generation time. The transformation to ECEF is done by first finding the conversion from TEME to True of Date. True Equator of Epoch this approach, the epoch of the TEME frame is held constant. Subsequent rotation matrices must therefore account for the change in precession and nutation from the epoch of the TEME. frame to the epoch of the transformation and no one really knows which one is used by NORAD when releasing TLEs [1] 1 , however, people generally believe the of date option is correct [1]. Regardless, when you compare the positions predicted by SGP4 to the positions stated in the ephemeris, you will need to make sure they are in the same coordinate frame! I recommend converting everything to the frame of the ephemeris data-set (which you should have identified in question 1). Using the software package or programming language and library of your choice (Excel, C, Python, MATLAB, STK, or GMAT) use SGP4 to obtain predicted locations for each of the time steps in the ephemeris data-set you were provided. For each time step, get the difference between the position predicted by SGP4, and the ephemeris data-set position. Plot the magnitude of these differences for each TLE that you downloaded. (15 pts) Question 4 In Figure 1, the errors between a blended set of TLEs and a Starlink ephemeris are shown. In particular, they are shown in the R , S,W frame, whose basis vectors are constructed as follows: R = r r , W = r v ?r v? , S = W R , (1) where r is the distance from the center of the earth to the orbit, r the position vector of the satellite, and v the velocity vector of the satellite. Both r and v are in the same frame as the ephemeris data-set; it is these new basis vectors that define our body centric frame. 1. Identify which basis vectors are radial, cross track, and along track. (3 pts) 2. For the first time step in your ephemeris data, calculate the 3 by 3 transformation matrix R that will map a relative point in the ephemeris coordinate frame to one in the R , S,W frame, and calculate RR? . (3 pts) Question 5 1Note the authors on this paper, and check out who they work for! 4 With the basis vectors constructed through the process above, we can move any relative position between an ephemeris time-step and, say, an SGP4 predicted position, into this body frame. Let rSGP 4 be a position predicted by SGP4 for a given point in time. Let rT be the ephemeris data-set position at that same point in time. Here, T implies truth, as we are considering this to be our ground truth. The difference vector dr = rSGP 4 - rT , when multiplied with our prior transformation matrix R will give a new position for rSGP 4, but in the R , S,W frame. Arithmetically: ?r = Rdr. (2) 1. For each time step, work out the relative error ?r between the ephemeris ground truth, and the SGP4 predicted position, for each TLE. If you do it right, you should get something similar to Figure 1, however growing in error with time. Remember to convert the outputs of the SGP4 propagation to the coordinate frame of the ephemeris dataset! (15 pts) 2. Which axis has the greatest error magnitude? Is it cross track, along track, or radial? (3 pts) 3. Give a rough estimate (no calculations necessary) in kilometers per day of the error growth rate of SGP4 according to your results. Use (and cite) a peer reviewed journal article to state if this performance is more/less accurate than would be expected from SGP4. (5 pts) 4. Suppose the error between SGP4 and the truth orbit was significantly greater than would be expected what could cause this, assuming you did everything else correctly? (3 pts) Question 6 3 Orbital Element Calculation and Perturbations From the previous section, you should now have a few sets of data the original ephemeris provided to you, and three SGP4 generated ephemerides. You will only need the original ephemeris and the data generated by the most up-to-date TLE. 1. For each time step in your ephemeris data-set plot each of the following (32 pts): semi-major axis a eccentricity ? inclination i right ascension of the ascending node ? first derivative of the right ascenscion of the ascending node ?? argument of perigee ? true anomaly ? Mean anomaly M 2. In theory, mean anomaly should progress smoothly, producing a saw-tooth like plot, going from 0 ? to 360? in the time it takes to do a single orbit. And in a very-close-to-circular orbit, like that which youve been given, true anomaly should do the same. Is this what your plots show? If they exhibit more complex behaviour, what could be the cause for that? (7 pts) 3. For each time step, sum together the true anomaly and argument of perigee values, take the floating point modulo 360 (fmod(x,360) in most languages) and plot them. What do you see now? Why does this now create such a perfect plot? (5 pts) Question 7 5 You should have observed in your plot of ? that it varies with time. This is due to the slightly non-spherical nature of the earth (if we exclude other factors, such as manoeuvres), and is dependent on proximity to the earth, and inclination of the satellite orbit. This property is often used by sunsynchronous remote sensing satellites to ensure that the sun is always behind them, and illuminating the earth for half their orbit, to ensure high quality images. With orbital state vectors that use a cartesian coordinate system, we use what are called linearized J2 perturbation equations. Based on the X, Y, Z position of the craft, a perturbation acceleration can be calculated with the following: aJ2 = -J2 3 2 Earth |r| 2 REarth |r| 2 ? ? ? ? ? ? ? ? ? ? ? ? X |r| 1 - 5 h Z |r| i2 Y |r| 1 - 5 h Z |r| i2 Z |r| 3 - 5 h Z |r| i2 ? ? ? ? ? ? ? ? ? ? ? ? (3) Note the strong dependence of the final value for aJ2 on the Z value of the spacecrafts position. In Figure 3, the Z axis is vertical, though the zonal harmonic pictured there is greater than that of J2. Regardless, this suggests a heavy dependence of J2 nodal precession on the position of the spacecraft above/below the equator. We can calculate the theoretical rate of change of ? using following equation ? =? - 3 2 J2 REarth l M? cosi, (4) where J2 is the zonal harmonic (see Figure 3 to see the difference between zonal, tesseral, and sectoral spherical harmonics) for Earth, REarth is. . . the radius of the Earth, l is the semi-latus rectum of the orbit, which you should have calculated as part of generating values for the previous question, M? is the mean motion of the satellite, and i the inclination of the satellites orbit. (a) (b) (c) Figure 3: Difference between zonal, tesseral, and sectoral harmonics. Note that the order of the zonal harmonic depicted here is higher than that of the J2 zonal, hence it has more ridges. 6 1. Based on your results for ?? in the previous question, obtain an average ?? . (2 pts) 2. For each time step, plot ?aJ2 ?. (5 pts) 3. Based on the equation for ?? that uses the J2 zonal harmonic, obtain a theoretical value for ?? , and show working out. Does it match the average from the previous step? (5 pts) 4. Plot the latitude of the satellite as a function of time, and compare that plot to the plot of ?aJ2 ?. What conclusions can you draw from this? (5 pts) Question 8 References [1] D. Vallado, P. Crawford, R. Hujsak, and T. Kelso, Revisiting spacetrack report# 3, in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2006, p. 6753. 4 Day of Year Conversion 1 void days2mdhms ( 2 int year , double days , 3 int & mon , int& day , int & hr , int & minute , double & sec 4 ) 5 { 6 int i , inttemp , dayofyr ; 7 double temp ; 8 // start indicies at 1 9 int lmonth [] = { 0, 31 , 28 , 31 , 30 , 31 , 30 , 31 , 31 , 30 , 31 , 30 , 31 }; 10 11 dayofyr = (int ) floor ( days ) ; 12 /* ----------------- find month and day of month ---------------- */ 13 if (( year % 4) == 0) 14 lmonth [2] = 29; 15 16 i = 1; 17 inttemp = 0; 18 while (( dayofyr > inttemp + lmonth [ i ]) && ( i < 12) ) 19 { 20 inttemp = inttemp + lmonth [ i ]; 21 i ++; 22 } 23 mon = i ; 24 day = dayofyr - inttemp ; 25 26 /* ----------------- find hours minutes and seconds ------------- */ 27 temp = ( days - dayofyr ) * 24.0; 28 hr = floor ( temp ) ; 29 temp = ( temp - hr ) * 60.0; 30 minute = floor ( temp ) ; 31 sec = ( temp - minute ) * 60.0; 32 } Listing 1: C++ day of year to mdhms conversion code, David Vallado 7