%astro - demo of solar system orbits
r_earth = 1;   r_mars = 1.523;   %Orbital radius in AU
P_earth = 1;   P_mars = r_mars^1.5;   %Orbital periods in Yrs using Kepler III

t = linspace(0, 20*P_mars, 1000 );  %Time array in years
theta_earth = 2*pi*rand + 2*pi*t/P_earth;   
theta_mars = 2*pi*rand + 2*pi*t/P_mars; %Array of plant orbit polar angles
x_earth = r_earth*cos(theta_earth); 
y_earth = r_earth*sin(theta_earth);   %x,y coordinates in AU
x_mars = r_mars*cos(theta_mars); 
y_mars = r_mars*sin(theta_mars);

plot( x_earth,y_earth,'b','linewidth',2); hold on;  %Plot the orbits
plot( x_mars,y_mars,'r','linewidth',2); fsize=18; set(gca,'fontsize',fsize);
xlabel('x /AU'); ylabel('y /AU'); title('Earth and Mars orbits'); axis equal; grid on; hold on;
pearth = plot( x_earth(1),y_earth(1),'b*','markersize',42 ); 
pmars = plot( x_mars(1),y_mars(1),'r*','markersize',42 );

stop = 1; n=1;   %Animate the orbits
while stop==0
    n = n+1; if n>length(t); n=1; end
    set( pearth, 'xdata',x_earth(n),'ydata',y_earth(n) );  %Updating the x,y coordinates
    set( pmars, 'xdata',x_mars(n),'ydata',y_mars(n) );
    pause(0.05);
    drawnow;
end

%Step through the orbital times and draw a line between the positions of
%earth and mars
for n=1:length(t)
    plot( [ x_earth(n),x_mars(n) ],  [ y_earth(n),y_mars(n) ], 'k-' );
end

figure;  %Plot orbit of mars from perspective of earth
plot( x_mars - x_earth,y_mars - y_earth ,'r','linewidth',2);
xlabel('x /AU'); ylabel('y /AU'); title('Mars orbit from Earth'); 
axis equal; grid on;

