1ASEN 3200 Orbit Mechanics / Attitude Dynamics & Control Lab O-3 Assigned 2/24/06, Due 3/7/06 The potential function and equations of motion for a satellite perturbed by Earth oblateness are given below. []ECIECI22ECI22222232XXYYZZJRZ v113,Er2r 2rˆˆ ˆˆˆˆr I J K XI YJ ZKXYZWhereX3R ZX1J512rrr∂Φ⎡⎤⎢⎥∂⎡⎤⎢⎥⎢⎥∂Φ⎢⎥=∇Φ =⎢⎥⎢⎥∂⎢⎥⎢⎥⎣⎦∂Φ⎢⎥⎢⎥∂⎣⎦⎡⎤⎛⎞µ⎛⎞Φ= + − = −Φ⎢⎥⎜⎟⎜⎟⎝⎠⎢⎥⎝⎠⎣⎦∂Φ ∂Φ ∂Φ=∆Φ= + + = + +∂∂ ∂⎡⎛⎞µ⎛⎞=− − −⎢⎜⎟⎜⎟⎝⎠⎢⎝⎠⎣  22232YYXXZ3R ZZ1J532rrr⎤⎥⎥⎦=⎡⎤⎛⎞µ⎛⎞=− − −⎢⎥⎜⎟⎜⎟⎝⎠⎢⎥⎝⎠⎣⎦  322J .00108263, 398600.4418km /s , R 6378.137km=µ= =. J2, µ, and R are parameters for the JGM-2 gravity field. E is the total energy and will be constant for a potential function dependent only on zonal harmonics. The initial conditions for an orbit are given as follows: P=6200.0 sec, e=0.001, i=45°, Ω(RAAN) = 10.0, ω = 45°, θ=315°. Use an epoch and initial time of 1 Jan 2000 00:00:00 and the J2000 coordinate system. Use a stop time that is three orbit periods later. Assume the orbit is perturbed only by Earth oblateness, J2. Perform the following study: 1. Write the appropriate Matlab code and numerically integrate the equations of motion given above in the ECI Cartesian frame for three revolutions. Compute and store the orbit elements and ECI Cartesian coordinates at 1 min intervals. Starting with the same2initial conditions, generate the ECI Cartesian coordinates for a two-body orbit for the same time interval. a. Plot the orbit elements, the energy deviation from its epoch value, and the deviation of the Z component of angular momentum from its epoch value for the perturbed orbit. Both the energy and hz should deviate from their epoch values by around 10-7 if your integration is done properly. Use the equation for energy given above. b. Generate a file of Cartesian position coordinate differences between the perturbed and two body orbits (rrpeturbed body−−2). c. Code the RIC transformation matrix and plot the RIC differences using the two-body orbit as the reference. 2. Generate the perturbed orbit using STK. You will need to choose HPOP (High Precision Orbit Propagator) from the satellite basic properties. After choosing HPOP, click on the force model and choose the JGM2 gravity model. Set the maximum degree to 2 and the maximum order to 0. This will select J2 as the perturbing force. Turn off all other perturbing forces (drag, solar pressure, lunar, and solar gravity). Generate the RIC plots between the perturbed orbit and the two-body orbit as the reference using STK. Your RIC plots should look identical to those from MATLAB. You may want to do the STK plots first so that you will know if your MATLAB results are correct. 3. Your report should contain the following plots a. Plots of the orbit elements (a, e, i, Ω, ω,Tp), and the energy and hz deviations for your perturbed orbit generated with MATLAB. Do not plot the true anomaly. b. RIC differences between your perturbed and two body orbits. c. RIC differences between STK’s perturbed and two-body orbit. 4. In addition to observations you wish to make, discuss the following topics in your report: a. Verify the differential equations of motion by generating the partial derivative, X∂Φ∂, by hand as requested in the pre-lab assignment. b. Describe the variations of the orbit elements and note the frequencies that are present, i. e., secular, short period, and long period. Can you detect any long period perturbations? Why or why not? Is the energy and hz constant? c. Discuss the RIC result. What is the direction of the dominant variation from the two-body orbit? Why does the major perturbation occur in this direction? d. Do the STK results differ from yours? If so, why?3 The RIC transformation for part 1 is obtained as follows (see also the handout on coordinate systems). Given the position and velocity vector of the reference orbit in the ECI frame, we have []xyzXYZxyzRICxy zrhˆˆˆ ˆˆˆˆˆˆˆˆrXiYjZk,rXiYjZk,hrxrthen,R ,C ,ICxR.rhIn matrix notation,ˆˆˆiiRRRRˆˆˆIIIIjTjˆˆˆCC CCkk=++ =++ = = = =⎡⎤ ⎡⎤⎡⎤⎡⎤⎢⎥ ⎢⎥⎢⎥⎢⎥==⎢⎥ ⎢⎥⎢⎥⎢⎥⎢⎥ ⎢⎥⎢⎥⎢⎥⎢⎥ ⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ ⎣⎦ Note that the elements of the RIC transformation matrix at each point are the components of the RIC unit vectors. Hence, []XYZRICRXITYCZ∆⎡⎤ ⎡ ⎤⎢⎥ ⎢ ⎥=∆⎢⎥ ⎢ ⎥⎢⎥ ⎢ ⎥∆⎣⎦ ⎣ ⎦ Pre-Lab assignment: Verify the equation for X by formingX∂Φ∂. PROGRAMING NOTES The ode45 integrator in matlab is a variable step integrator. This means that it chooses its own step size based on default tolerances built into the program. These tolerances are very loose and as a result ode45 gives very poor results when integrating the orbit equations of motion. To solve this problem you can reset the tolerances by replacing the [t,x] statement with the following. tol = 1e-10; options = odeset('RelTol',tol,'AbsTol',tol); [t,x] = ode45('YourDerivativeFunction', time, x0, options) the time array can be defined as time=(0:60:18600) Where 0 is the initial time, 60s is the time step at which the ephemeris will be returned and 18600 sec is the final