Orbits with Kepler's Equation and Gravity Term J2

The two sets of plots below compare calculated with Kepler's equation with orbits calculated with gravity term J2.

Figure 1: Output from Kepler's equation
Figure 2: Orbits calculated from J2
1. % demoOrbit
2. % This demo compares the orbits due to Kepler's equations
3. % that assume that the earth is a point mass, and
4. % the orbit of a satellite including the effect of the
5. % gravity term J2 (flattening of Earth at the poles)
6. %
7. % Initial condition
8. % vx vy vz x y z
9. X0 = [0 -80007200]';
10. T = 0:10:16000;
11. lsode_options('absolute',1e-10);
12. lsode_options('relative',1e-10);
13. Y0 = lsode('satKepler', X0, T, 0);
14. R0 = Y0(:,4:6)';
15. V0 = Y0(:,1:3)';
16. h  = cross(R0,V0);
17. Adot = 1/2*sqrt(h(1,:).*h(1,:) + h(2,:).*h(2,:) + h(3,:).*h(3,:));
18. YJ2  = lsode('satJ2', X0, T, 0);
19. RJ2  = YJ2(:,4:6)';
20. VJ2  = YJ2(:,1:3)';
21. h    = cross(RJ2,VJ2);
22. AdotJ2 = 1/2*sqrt(h(1,:).*h(1,:) + h(2,:).*h(2,:) + h(3,:).*h(3,:));
23. x    = Y0(:,4);
24. y    = Y0(:,5);
25. z    = Y0(:,6);
26. vx   = Y0(:,1);
27. vy   = Y0(:,2);
28. vz   = Y0(:,3);
29. r    = sqrt(x.*x + y.*y + z.*z);
30. v    = sqrt(vx.*vx + vy.*vy + vz.*vz);
31. xJ2  = YJ2(:,4);
32. yJ2  = YJ2(:,5);
33. zJ2  = YJ2(:,6);
34. vxJ2 = YJ2(:,1);
35. vyJ2 = YJ2(:,2);
36. vzJ2 = YJ2(:,3);
37. rJ2  = sqrt(xJ2.*xJ2 + yJ2.*yJ2 + zJ2.*zJ2);
38. vJ2  = sqrt(vxJ2.*vxJ2 + vyJ2.*vyJ2 + vzJ2.*vzJ2);
39. erase
40. figure(1)
41. window('221')
42. plot(T,r,'blue',T,rJ2,'red');
43. ylabel('r, rJ2 [km]')
44. window('222')
45. plot(T,v,'blue',T,vJ2,'red');
46. ylabel('v [km/s]')
47. window('223')
50. window('224')
51. plot(y,z,'blue',yJ2,zJ2,'red','grid');
52. ylabel('z, zJ2 [km]')
53. xlabel('y, yJ2 [km]')
54. figure(2)
55. window('221')
56. plot(T,rJ2-r,'red');
57. ylabel('rJ2-rKepler [km]')
58. window('222')
59. plot(T,vJ2-v,'red');
60. ylabel('vJ2-vKepler [km/s]')
61. window('223')