MLDesign Technologies, Inc.

Mainnavigation

Subnavigation

BORDER

 
 
 

Pagecontent

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')
  48. plot(T,Adot,'blue',T,AdotJ2,'red')
  49. ylabel('Adot, AdotJ2')
  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')
  62. plot(T,AdotJ2-Adot,'red')
  63. ylabel('AdotJ2-AdotKepler')
  64. window('224')
  65. plot(yJ2-y,zJ2-z,'red','grid');
  66. ylabel('zJ2-z [km]')
  67. xlabel('yJ2-y [km]')