% Juno - program to compute the trajectory of the Juno probe to Jupiter
% The probe is launched from a low earth orbit into a two year orbit around
% the sun.  A guidance boost is applied approximately one year into the
% flight.  The probe returns to earth after two years and the slingshot
% effect accelerates the probe toward Jupiter.

clear
G=6.7e-11;            % gravitational constant
mSun = 1.9e30;        % mass of the sun
rEarthOrbit = 150e9;  % radius of earth orbit around sun
mEarth = 5.9742e24;   % mass of the earth
rEarth = 6.378e6;     % radius of earth

% parameters for earth orbit
xE(1)=rEarthOrbit; 
vxE(1)=0; 
yE(1)=0; 
vyE(1)=28500;         % velocity of earth 

% launch parameters for probe 
% from a low earth orbit with altitude = 400 km, velocity = 7000 m/s
xR(1)=rEarthOrbit + rEarth + 400000; 
yR(1)=0;
LV=4900;  % the probe launch velocities were determined by trial and error 
vxR(1)=3000;   % to offset earth's gravity after launch      
vyR(1)=28500 + 7000 + LV; % total initial probe velocity 
                      % = earth orbit vel. + probe orbit vel. + launch vel.

% midcourse  boost parameters, again, lots of trial and error
TBstart = 345*60*60*24;  % apply boost 345 days into flight
TB = 5*60; % apply boost for 5 minutes  
BOOSTy = 3.1017; % debug - 3 late, 3.25 early, 3.1  slingshot

dREmin = 0;   % the rocket / earth point of closest approach

t(1) = 0; % start time
dt = 60;  % dt is too big for accuracy, to keep run time reasonable
          % if smaller the boost numbers will have to be changed
n = 3.1*365*60*60*24/dt;  % fly for 3.1 years

for i=1:n
    t(i+1) = t(i) + dt;          
    xE(i+1) = xE(i) + vxE(i)*dt;  % update earth position
    yE(i+1) = yE(i) + vyE(i)*dt;  
    
    dES = sqrt(xE(i)^2+yE(i)^2);  % calculate earth accel    
    AE = G*mSun/dES^2;            
 
    vxE(i+1)= vxE(i) - AE*(xE(i))/dES*dt;  % update earth velocity
    vyE(i+1)= vyE(i) - AE*(yE(i))/dES*dt;  

    xR(i+1) = xR(i) + vxR(i)*dt;  % update rocket position
    yR(i+1) = yR(i) + vyR(i)*dt;  
    
    dRS = sqrt(xR(i)^2+yR(i)^2);  % calculate rocket accel due to sun   
    ARS = G*mSun/dRS^2;            
 
    dRE = sqrt((xR(i)-xE(i))^2+(yR(i)-yE(i))^2); % calculate rocket accel  
    ARE = G*mEarth/dRE^2;                        % due to earth  
    
    % update rocket velocity
    vxR(i+1)= vxR(i) - (ARS*xR(i)/dRS + ARE*(xR(i) - xE(i))/dRE)*dt;  
    vyR(i+1)= vyR(i) - (ARS*yR(i)/dRS + ARE*(yR(i) - yE(i))/dRE)*dt;  
    
    % apply midcourse guidance boost
    if (t(i) >= TBstart && t(i) < TBstart + TB)
      vyR(i+1) = vyR(i+1) + BOOSTy*dt;     % add accel due to boost
      dREmin = dRE;      % initialize point of closest approach on return
   end 

   if dRE < dREmin     % record point of closest approach for tweaking
       dREmin = dRE;   % boost
       dREmini = i;
   end

end

II = 1:i+1;                  % whole trajectory
%II = RERi-1000:RERi+1000;   % flyby trajectory
figure(1)
plot(xR(II),yR(II),'r'); hold on; plot(xE(II),yE(II),'b'); hold off 
axis equal
title('Juno trajectory, earth blue, probe red')
xlabel('x')
ylabel('y')
zlabel('years')

figure(2)                 % x,y,t graph
plot3(xR(II),yR(II),t(II)/(365*24*60*60), 'red')
hold on
plot3(xE(II),yE(II),t(II)/(365*24*60*60), 'blue')
grid
title({'Juno trajectory';'earth - blue, probe - red'})
xlabel('x,')
ylabel('y')
zlabel('years')
hold off


