% Program to compute 2D orbits.

clear
G = 6.7000e-11;  % gravitational constant
Mearth = 5.9e24;  % mass of earth, kg
Rearth = 6.3e6;  % radius of earth, m

% orbit initial conditions 
% moon: p = (0, Rearth+380000000), v = (1000, 0) period = 30 days
% space station: (0, Rearth+402000), (7706, 0), 93 minutes 
% ellipse: (0, Rearth+402000), (10000, 0), 600 minutes
% hyperbola: (40000000, 12000000), (-10000, 0), 110 minutes 

px(1)= 0;              % initial position (for space station)
py(1)= Rearth+402000;    
vx(1)=7706 ;           % initial velocity (for space station) 
vy(1)= 0; 

t(1)=0;                % start time
dt=.1;

for i=1:93*60/dt;   % run for one period
    t(i+1)=t(i) + dt;
    
    px(i+1)=px(i) + vx(i)*dt;  % update position
    py(i+1)=py(i) + vy(i)*dt;  
    
    R=sqrt(px(i)^2 + py(i)^2);  % calculate acceleration
    A= -G*Mearth/R^2;           
    Ax= A*px(i)/R;              % resolve into x and y components
    Ay= A*py(i)/R;              
    
    vx(i+1)= vx(i) + Ax*dt;  % update velocity
    vy(i+1)= vy(i) + Ay*dt;  
end

plot(px, py)         % plot satellite trajectory
hold on              % draw circle for earth
for j = 1:361   
    theta = j*(2*pi/360);   
    cx(j)= cos(theta);
    cy(j)= sin(theta);
end
plot(cx*Rearth,cy*Rearth)  % draw earth
xlabel('x')
ylabel('y')
title('Space station orbit')
%title('Elliptical orbit')
%title('Hyperbolic trajectory')
axis('equal')
hold off

    
    
