% Program to shoot a rocket to the moon.  We launch the rocket straight
% up from the north pole, i.e. along the y axis.  First we determine by
% trial and error the necessary launch velocity to just get the rocket up
% to the moon's orbit.  Then we position the moon so that it reaches the
% point above the rocket launch point just as the rocket reaches the moon's
% orbit.  Then we tweak the velocity to get a hit.

G=6.7e-11;                % gravitational constant, m3/kg/s2
mEarth = 5.9742e24;       % mass of the earth, kg
rEarth = 6.378e6;         % radius of the earth, m

rMoon = 1.737e6;          % radius of moon, m
mMoon = 7.35e22;          % radius of moon. kg
rMoonOrbit = 385000000;   % radius of moon's orbit, m

period = 27.3*24*60*60;       % one sidereal month in seconds
vM=pi*2*(rMoonOrbit)/period;  % calcululate moon velocity

% determine rocket velocity to just reach moon orbit by trial and error
vM0 = 11170;     % velocity to just get rocket to moon orbit

% calculate initial position of moon so it will be above rocket launch 
% point when rocket reaches moon orbit, we could do the calculations 
% offline if we wanted
DAYS = 3.9;      % observe time it takes for rocket to reach moon orbit
theta = (DAYS*24*60*60 / period) * 2 * pi; % moon's position off vertical
alpha = pi/2 - theta;   % moon's position off horizontal

vM0 = 11170 - 55 ; % tweak rocket velocity to hit moon, by trial and error
%VEL = 11170 - 55.17; % tweak to just miss for Apollo, by trial and error

% initial parameters for moon orbit, theta radians off y axis
px(1) = rMoonOrbit*cos(alpha);     % initial position 
py(1) = rMoonOrbit*sin(alpha);   
vx(1) = vM*cos(alpha+pi/2);   % rotate velocity vector(0,vM)alpha radians
vy(1) = vM*sin(alpha+pi/2);      

% initial parameters for rocket
xR(1) = 0;        % initial x position
yR(1) = rEarth;   % launch from earth's surface
vxR(1) = 0;       % initial x velocity
vyR(1) = vM0;   

dt=1;
t(1) = 0;                 % start time
n = 4.4*24*60*60 / dt;    % fly for 3.9 days

for i=1:n
    t(i+1) = t(i) + dt;          
    px(i+1) = px(i) + vx(i)*dt;  % update moon position                     
    py(i+1) = py(i) + vy(i)*dt;  
    
    Rmr = sqrt(px(i)^2+py(i)^2);  % calculate moon acceleration    
    Amr = -G*mEarth/Rmr^2;            
 
    vx(i+1)= vx(i) + Amr*(px(i))/Rmr*dt;  % update moon velocity
    vy(i+1)= vy(i) + Amr*(py(i))/Rmr*dt;  

    xR(i+1) = xR(i) + vxR(i)*dt;  % update rocket position
    yR(i+1) = yR(i) + vyR(i)*dt;  
    
    Rr = sqrt(xR(i)^2+yR(i)^2);  % calculate rocket accel due to earth    
    Are = -G*mEarth/Rr^2; 
    
    Rmr = sqrt((xR(i)-px(i))^2+(yR(i)-py(i))^2); % calculate rocket accel 
    Amr = -G*mMoon/Rmr^2;          % due to moon  
 
    % update rocket velocity 
    vxR(i+1)= vxR(i) + Amr*(xR(i) - px(i))/Rmr*dt + Are*xR(i)/Rr*dt; 
    vyR(i+1)= vyR(i) + Amr*(yR(i) - py(i))/Rmr*dt + Are*yR(i)/Rr*dt;   

    if (Rmr < rMoon)   % has rocket hit moon?
        'hit'          % yes
        break 
    end; 
end

figure(1)
II=[i+1-40000:i+1];  % graph last 40000 seconds of rocket/moon trajectories
plot(px(II),py(II))  % plot moon trajectory
hold on
plot(xR(II),yR(II))  % plot rocket trajectory

for j=1:101          % draw the moon centered at (px(i),py(i))
   ang=2*pi*j/100;
   xG(j)= rMoon*cos(ang);
   yG(j)= rMoon*sin(ang);
end
plot(xG + px(i),yG + py(i))  

xlabel('x')
ylabel('y')
title('Shoot the moon')
%title('Slingshot')
hold off

figure(2)               % plot velocity profiles for slingshot
plot(t(II),sqrt(vxR(II).^2+vyR(II).^2),'blue')
hold on
plot(t(II),vxR(II),'red')
plot(t(II),vyR(II),'green')
xlabel('t')
ylabel('v in m/s')
title('vR - blue, vxR - red, vyR - green')
hold off



