Background The Gravity Recovery And Climate Experiment (GRACE) mission was launched by NASA in 2002. The primary objective of GRACE was to measure variations in Earth’s gravity field to track mass changes in the hydrosphere, cryosphere,solid Earth, and oceans. In contrast to single satellite approaches with one dedicated sensor, GRACE used a constellation of two satellites, orbiting one behind the other, featuring a suite of measurement systems. GRACE mission ended in July 2017. The GRACE FollowOn (GRACE-FO) mission, which launched May 22, 2018, will continue the work of GRACE mission. GRACE and GRACE-FO are very similar satellite missions. Each satellite mission consists of two identical satellites flying in the same near-polar orbit at about 500 km altitude, one following the other at a distance of about 220 km (see Figure 1 below). Both satellites are equipped with (1) GPS to measure their position and (2) accelerometers (placed at the center of mass of the satellite) to measure nongravitational forces like air drag acting on the satellites. The principal observation of the GRACE and GRACE-FO missions is the ultra-precise measurement of inter-satellite distance between the two satellites. Distance between the two satellites changes because of the gravitational pull of the masses beneath the satellites. Therefore, tracking the change in the inter-satellite distance enables us to measure mass changes that happen on or near the Earth’s surface. Surface mass changes are caused by various processes like the Greenland icesheet melting due to global warming and climate change or groundwater depletion in California’s Central Valley due to overpumping of groundwater used for irrigation. GRACE-FO Orbit Data from NASA Jet Propulsion Laboratory (JPL) The files “GNV1B_2024-02-22_C_04.txt” and “GNV1B_2024-02-22_D_04.txt” provide the orbits of GRACE-FO satellite #1 and #2 (also known as GRACE-FO satellite C and D) on 22 February 2024, respectively, released by NASA Jet Propulsion Laboratory (JPL). The files include the position and velocity vector (as well as their error information) of each satellite every second and the time (epoch) of position/velocity vectors. The time is given as seconds passed since 01-Jan-2000 11:59:47 UTC. Position and velocity vectors are given in the Earth-Centered Earth-Fixed (ECEF) geocentric coordinate system: The Z axis is the line between the North and South Poles (with positive values increasing northward), the X axis is in the plane of the equator, passing through the origin and prime meridian (Greenwich), and the Y axis in the plane of the equator, extending from 90°W longitude to 90°E longitude (positive). The MATLAB function “ReadGFO_Orbit.m” provided to you reads the information given in the GRACE-FO orbit files and outputs the time (in decimal year) as well as position and velocity vector of the GRACE-FO satellites #1 and #2 (also known as satellites C and D). You will use this function to extract the information required for the programming tasks described below from the GRACE-FO files. ASEN 1320 Aerospace Computing and Engineering Applications Spring 2025 Figure 1. GRACE and GRACE-FO measurement is implemented by two identical satellites orbiting one behind the other in a near-polar orbit plane at ~500 km altitude. The along-track separation is kept within a range of 220 ± 50 km. The satellites experience positive and negative, gravitationally induced, along-track accelerations due to the varying mass distribution below them. Each satellite will experience the effects of the local mass at slightly different times causing a differential acceleration. The differential acceleration, in turn, causes distance (range) variations and velocity differences Δv that are proportional to the mass attraction related to the gravitational potential, U. The relative distance between the satellites is measured with micron-level precision by a high accuracy dual frequency K-band inter-satellite ranging system. An accurate three-axis accelerometer measures the effects of all non-gravitational forces acting on each satellite, including atmospheric drag, direct and Earth-reflected solar radiation pressure, and thrusting. A GPS receiver on each satellite provides position and time synchronization. The satellites overfly the entire Earth surface within approximately 30 days, allowing monthly estimates of a global gravity model with a surface spatial resolution of about 300 km with an accuracy of 2 cm. [Reference: Tapley et al., 2019]. ASEN 1320 Aerospace Computing and Engineering Applications Spring 2025 Programming Tasks Write six user-defined MATLAB functions and a main program to perform the programing tasks described below. All the user-defined functions should be called inside the main function to perform the required actions. [5 points for calling the functions in the main program] 1. Review the function called ReadGFO_Orbit (provided to you) with the following input and output parameters: Input  A string variable called InFilename, indicating the name of the GRACE-FO orbit text file (for example, “GNV1B_2024-02-22_C_04.txt”). Output  A column vector called t_dy, indicating the time (epoch) of satellite position/velocity measurements converted to decimal year.  A matrix (with 3 columns) called r, indicating the position [x y z] of GRACE-FO satellite. The number of rows in this matrix is equal to the number of position measurements of GRACE-FO satellites given in text file (which is equal to 86400: 1 position per second).  A matrix (with 3 columns) called v, indicating the velocity [Vx Vy Vz] of GRACE-FO satellite. The number of rows in this matrix is equal to the number of velocity measurements of GRACE-FO satellites given in the text file (which is equal to 86400). This function will read the data from the text files “GNV1B_2024-02-22_C_04.txt” and “GNV1B_2024-02-22_D_04.txt” and extract time (and will convert it to decimal year) and position and velocity vector of the GRACE-FO satellite #1 and #2. Therefore, this function should be called twice in the main program. The MATLAB function “ReadGFO_Orbit.m” (provided to you) performs this programming task for you. Please go through this function carefully to make sure that you fully understand the MATLAB built-in functions used in this user-defined function. 2. Write a function called GFO_Range with the following input and output parameters: [10 points] Input  A matrix (with 3 columns) called r1, indicating the position [x1 y1 z1] of the GRACE-FO satellite #1.  A matrix (with 3 columns) called r2, indicating the position [x2 y2 z2] of the GRACE-FO satellite #2. Output  A column vector called Rho, indicating the inter-satellite range (distance) between the two GRACE-FO satellites. This function will compute the inter-satellite range (distance) 𝜌𝜌 using the following equation: 𝜌𝜌 = �(𝑥𝑥2 − 𝑥𝑥1)2 + (𝑦𝑦2 − 𝑦𝑦1)2 + (𝑧𝑧2 − 𝑧𝑧1)2 3. Write a function called GFO_RangeRate with the following input and output parameters: [15 ASEN 1320 Aerospace Computing and Engineering Applications Spring 2025 points] Input  A matrix (with 3 columns) called r1, indicating the position [x1 y1 z1] of the GRACE-FO satellite #1.  A matrix (with 3 columns) called r2, indicating the position [x2 y2 z2] of the GRACE-FO satellite #2.  A matrix (with 3 columns) called v1, indicating the velocity [Vx1 Vy1 Vz1] of the GRACE-FO satellite #1.  A matrix (with 3 columns) called v2, indicating the velocity [Vx2 Vy2 Vz2] of the GRACE-FO satellite #2.  Output  A column vector called Rho_dot, indicating the inter-satellite range-rate. This function will compute the inter-satellite range-rate 𝜌𝜌̇ using the following equations: 𝜌𝜌̇ =< 𝚫𝚫𝒗𝒗 ⋅ 𝒆𝒆𝟏𝟏𝟏𝟏 𝐋𝐋𝐋𝐋𝐋𝐋 > where < ⋅ > indicates the dot product and 𝚫𝚫𝒗𝒗 and 𝒆𝒆𝟏𝟏𝟏𝟏 𝐋𝐋𝐋𝐋𝐋𝐋 are defined as 𝚫𝚫𝒗𝒗 = 𝒗𝒗𝟐𝟐 − 𝒗𝒗𝟏𝟏 = [𝑉𝑉𝑥𝑥2 − 𝑉𝑉𝑥𝑥1 𝑉𝑉𝑦𝑦2 − 𝑉𝑉𝑦𝑦1 𝑉𝑉𝑧𝑧2 − 𝑉𝑉𝑧𝑧1] 𝒆𝒆𝟏𝟏𝟏𝟏 𝐋𝐋𝐋𝐋𝐋𝐋 = 𝚫𝚫𝒓𝒓 𝜌𝜌 = 𝒓𝒓𝟐𝟐 − 𝒓𝒓𝟏𝟏 𝜌𝜌 = [ 𝑥𝑥2 − 𝑥𝑥1 𝜌𝜌 𝑦𝑦2 − 𝑦𝑦1 𝜌𝜌 𝑧𝑧2 − 𝑧𝑧1 𝜌𝜌 ] 𝚫𝚫𝒓𝒓 = 𝒓𝒓𝟐𝟐 − 𝒓𝒓𝟏𝟏 = [𝑥𝑥2 − 𝑥𝑥1 𝑦𝑦2 − 𝑦𝑦1 𝑧𝑧2 − 𝑧𝑧1] 𝚫𝚫𝒗𝒗 is the relative velocity vector between the two satellites and 𝒆𝒆𝟏𝟏𝟏𝟏 𝐋𝐋𝐋𝐋𝐋𝐋 is the unit vector pointing from GRACE-FO satellite #1 to satellite #2. 4. Write a function called GFO_NUmDiff with the following input and output parameters: [15 points] Input  A scalar variable called dt, indicating the sampling rate (in seconds) of position and velocity vectors given in the GRACE-FO orbit files (which is equal to 1 second).  A column vector called Rho, indicating the inter-satellite range (distance) between the two GRACE-FO satellites.  A column vector called Rho_dot, indicating the inter-satellite range-rate. Output  A column vector called Rho_dot_ND, indicating the inter-satellite range-rate computed from inter-satellite range 𝜌𝜌 using numerical differentiation (see 𝜌𝜌̇ 𝑁𝑁𝑁𝑁 below).  A column vector called Rho_dot_diff, indicating the difference between Rho_dot_ND and Rho_dot (= 𝜌𝜌̇ 𝑁𝑁𝑁𝑁 − 𝜌𝜌̇). This function will compute the inter-satellite range-rate from inter-satellite range 𝜌𝜌 using numerical differentiation as follows (𝑁𝑁 is equal to the number of epochs (86400) in the GRACE-FO orbit files): ASEN 1320 Aerospace Computing and Engineering Applications Spring 2025 𝜌𝜌̇𝑁𝑁𝑁𝑁(𝑡𝑡𝑖𝑖) = 𝜌𝜌̇𝑁𝑁𝑁𝑁(𝑡𝑡1) = 𝜌𝜌(𝑡𝑡2) − 𝜌𝜌(𝑡𝑡1) 𝑑𝑑𝑑𝑑 𝜌𝜌(𝑡𝑡𝑖𝑖+1) − 𝜌𝜌(𝑡𝑡𝑖𝑖−1) 2𝑑𝑑𝑑𝑑 ; 𝑖𝑖 = 2, 3, 4, … , 𝑁𝑁 − 1 𝜌𝜌̇𝑁𝑁𝑁𝑁(𝑡𝑡𝑁𝑁) = 𝜌𝜌(𝑡𝑡𝑁𝑁) − 𝜌𝜌(𝑡𝑡𝑁𝑁−1) 𝑑𝑑𝑑𝑑 This function will also compute Rho_dot_diff which is equal to 𝜌𝜌̇𝑁𝑁𝑁𝑁 − 𝜌𝜌̇. 5. Write a function called SatVisibility with the following input and output parameters: [15 points] Input  A 1×3 row vector called rG, indicating the ECEF position [xG yG zG] of the Satellite Laser Ranging (SLR) station in Maryland, U.S.A. rG = [1130714.219 -4831369.903 3994085.962].  A matrix (with 3 columns) called r1, indicating the position [x1 y1 z1] of the GRACE-FO satellite #1. Output  A column vector called el, indicating the elevation angle (in degrees) of the GRACE-FO satellite #1 from the SLR station in Maryland.  A column vector called index_vis, indicating the index of elevation angles greater than 10 degrees in the column vector el. Such elevation angles indicate the epochs of GRACE-FO satellite #1 visible from the SLR station. The SLR technique measures the distance between the station and satellites using the laser pulses when the satellite is visible from the station. This function will first compute the elevation angle of the GRACE-FO satellite #1 with respect to the SLR station in Maryland as follows: 𝑒𝑒𝑒𝑒 = arcsin � 𝑧𝑧1 − 𝑧𝑧𝐺𝐺 𝑑𝑑 � 𝑑𝑑 = �(𝑥𝑥1 − 𝑥𝑥𝐺𝐺)2 + (𝑦𝑦1 − 𝑦𝑦𝐺𝐺)2 + (𝑧𝑧1 − 𝑧𝑧𝐺𝐺)2 where [𝑥𝑥1 𝑦𝑦1 𝑧𝑧1] and [𝑥𝑥𝐺𝐺 𝑦𝑦𝐺𝐺 𝑧𝑧𝐺𝐺] indicate the position of GRACE-FO satellite #1 and SLR station, respectively. The function will also return the index of elevation angle values index_vis in the column vector el that are greater than 10 degrees. 6. Write a function called writeGFO_CSV with the following input parameters: [10 points] Input  A string variable called OutFilename, indicating the name of the csv file with 3 columns which will be used to store the values of column vectors t_dy, Rho, and Rho_dot.  A column vector called t_dy, indicating the time (epoch) of satellite position/velocity measurements in decimal year.  A column vector called Rho, indicating the inter-satellite range (distance).  A column vector called Rho_dot, indicating the inter-satellite range-rate. Output: No output arguments. This function will write the values in column vectors t_dy, Rho, and Rho_dot into 3 columns of the csv file OutFilename. Use "time [dec. year]", "range [m]" and "range-rate [m/s]" as the header for each column. Call this function in the main program and use “GFO_file.csv” as the name of the output csv file. ASEN 1320 Aerospace Computing and Engineering Applications Spring 2025 7. A) In the MATLAB main program, write a code that produces a figure with 3 subplots displaying the time series of Rho [in km], Rho_dot [in m/s] and el [in degrees] with respect to time (given as hours on 22 February 2024). The elevation angles greater than 10 degrees should be marked with a different color. [15 points] B) In the MATLAB main program, write a code that produces a figure showing the ground track of GRACE-FO satellite #1 overlaid on a spatial map of the world. Portion of the ground track visible from the SLR station in Maryland with elevation angle greater than 10 degrees should be marked with a different color. Location of the SLR station should be shown with a black triangle. [15 points] To display the ground track of a satellite, the ECEF Cartesian coordinates [𝑥𝑥 𝑦𝑦 𝑧𝑧] of the satellite should be converted into the curvilinear coordinates [𝜑𝜑 𝜆𝜆 ℎ], where 𝜑𝜑 is latitude, 𝜆𝜆 is longitude and ℎ is the altitude from the Earth’s spherical surface. TIP #1: Use the MATLAB function “ecef2lla” to convert from Cartesian to curvilinear coordinates. TIP #2: Use the MATLAB function “geoplot” with latitude and longitude of the GRACE-FO satellite #1 to create this plot.