kepler1.m

航天工程工具箱

%KEPLER1  Demo of true, excentric and mean anomaly.
%   Demonstration of the relationship between the excentric and true anomaly.
%   The orbit of the object is correct in postion but incorrect
%   in velocity! The excentric anomaly is the angle in the circle to a
%   "virtual object" having the same position along the semi major axis of the
%   ellipse as the object following that orbit, whilst the mean anomaly is the mean
%   angle the orbit has as it would describe a circle with the same orbital period as
%   the ellipse.

% Copyright (c) 2001-10-27, B. Rasmus Anthin.

a=.5;
e=default(input('Enter eccentricity [0.8]: '),.8);
c=a*e;
b=sqrt(a^2-c^2);
N=default(input('Enter number of points [250]: '),250);
loop=default(input('Enter number of revolts [3]: '),3);
M=linspace(0,2*pi,N);
Ef=inline('M-E+e*sin(E)','E','e','M');
clear E
warning off
E(1)=fzero(Ef,M(1),2/N,0,e,M(1));
for i=2:N
   E(i)=fzero(Ef,E(i-1),2/N,0,e,M(i));
end
warning on
theta=unwrap(2*atan(sqrt((1+e)/(1-e))*tan(E/2)));

hold on
plot(a*cos(E),a*sin(E))
plot(a*cos(E),b*sin(E))
plot(-c,0,'rx',c,0,'rx')
plot([-a a],[0 0],'g:')
h=[];
axis(.5*[-1 1 -1 1])
axis square
for j=1:loop
   for i=1:length(M)
      delete(h)
      r=a*(1-e^2)./(1+e*cos(theta(i)));
      x=r*cos(theta(i));y=r*sin(theta(i));
      h=plot(c+[0 x],[0 y],'k',...                   %true anomaly
         c+x,y,'ko',...
         c+.05*cos(theta(1:i)),.05*sin(theta(1:i)),'b',...
         [0 a*cos(E(i))],[0 a*sin(E(i))],'k',...     %eccentric anomaly
         a*cos(E(i))*[1 1],[a*sin(E(i)) 0],'k:',...
         .05*cos(E(1:i)),.05*sin(E(1:i)),'b',...
         [0 a*cos(M(i))],[0 a*sin(M(i))],'r',...     %mean anomaly
         .08*cos(M(1:i)),.08*sin(M(1:i)),'m');
      drawnow
   end
end
figure,plot(M,theta,M,E),xlabel M,legend('\theta','E')