% Sedivy ECE 6390 Project Orbital Path 08JUL13
% Earth Orbit- Venus Orbit- Sat Interplanetary Transition Orbit
clear all
GMs = 1.3274935e11;   % Kepler's Constant for the Sun (km^3/s^2)
Max_Iteration = 6e4;
dt = 600;             % time increment (seconds)

% Plot for Earth Orbit
% Earth Initial Conditions
rearth = 152098232;         % initial radius @ aphelion (km)
theta_earth = 0;            % inital angle (by definition) (radians)
v_r_earth = 0;              % initial r-component of velocity (zero at aphelion)(km/s)
v_theta_earth = 29.295078;  % initial theta-component of velocity @ aphelion (km/s)

% Discrete Parameters and initial values for Earth
dr_earth = v_r_earth*dt;                  % km
dtheta_earth = v_theta_earth/rearth*dt;  % radians
rearth_0 = rearth;
theta_earth_0 = theta_earth;

% Set up arrays for r and theta
all_rearth = rearth;
all_theta_earth = theta_earth;
     
while abs(theta_earth - theta_earth_0) <= 2*pi && length(all_rearth) < Max_Iteration
    % Compute new state information
    rearth_new = rearth + dr_earth;
    theta_earth_new = theta_earth + dtheta_earth;
    dr_earth_new = dr_earth + ((rearth + 0.5*dr_earth)*dtheta_earth^2 - GMs/(rearth^2)*dt^2);
    dtheta_earth_new = dtheta_earth - (2*dr_earth*dtheta_earth/(rearth + 0.5*dr_earth));
    
    % Store information for r and theta
    all_rearth = [all_rearth, rearth_new];
    all_theta_earth = [all_theta_earth, theta_earth_new];

    % Change new state information to current state information for next iteration
    rearth = rearth_new;
    dr_earth = dr_earth_new;
    theta_earth = theta_earth_new;
    dtheta_earth = dtheta_earth_new;
end;

% Plot for Venus Orbit
% Venus Initial Conditions
rvenus = 107477000;         % initial radius @ perihelion (km)
theta_venus = 0;            % inital angle (by definition) (radians)
v_r_venus = 0;              % initial r-component of velocity (zero at perihelion)(km/s)
v_theta_venus = 35.263114;  % initial theta-component of velocity @ perihelion (km/s)

% Discrete Parameters and initial values for venus
dr_venus = v_r_venus*dt;                  % km
dtheta_venus = v_theta_venus/rvenus*dt;  % radians
rvenus_0 = rvenus;
theta_venus_0 = theta_venus;

% Set up arrays for r and theta
all_rvenus = rvenus;
all_theta_venus = theta_venus;
     
while abs(theta_venus - theta_venus_0) <= 2*pi && length(all_rvenus) < Max_Iteration
    % Compute new state information
    rvenus_new = rvenus + dr_venus;
    theta_venus_new = theta_venus + dtheta_venus;
    dr_venus_new = dr_venus + ((rvenus + 0.5*dr_venus)*dtheta_venus^2 - GMs/(rvenus^2)*dt^2);
    dtheta_venus_new = dtheta_venus - (2*dr_venus*dtheta_venus/(rvenus + 0.5*dr_venus));
    
    % Store information for r and theta
    all_rvenus = [all_rvenus, rvenus_new];
    all_theta_venus = [all_theta_venus, theta_venus_new];

    % Change new state information to current state information for next iteration
    rvenus = rvenus_new;
    dr_venus = dr_venus_new;
    theta_venus = theta_venus_new;
    dtheta_venus = dtheta_venus_new;
end;

% Plot for Satellite Transition Orbit
% Satellite Initial Conditions
daysoffset = -87.0;
valueoffset = floor(daysoffset*24*6); % put into units of dt, ensure it is an integer
if valueoffset < 0
    valueoffset = length(all_theta_earth) + valueoffset;
elseif valueoffset > length(all_theta_earth)
    valueoffset = valueoffset - length(all_theta_earth);
end;
rsat = all_rearth(valueoffset);            % initial radius @ launch (km)
theta_sat = all_theta_earth(valueoffset);  % inital angle @ launch (radians)
v_r_sat = 0;              % initial r-component of velocity (zero at aphelion)(km/s)
v_theta_sat = 27.270416;  % initial theta-component of velocity @ aphelion (km/s)

% Discrete Parameters and initial values for sat
dr_sat = v_r_sat*dt;                  % km
dtheta_sat = v_theta_sat/rsat*dt;  % radians
rsat_0 = rsat;
theta_sat_0 = theta_sat;

% Set up arrays for r and theta
all_rsat = rsat;
all_theta_sat = theta_sat;
     
while abs(theta_sat - theta_sat_0) <= 2*pi && length(all_rsat) < Max_Iteration
    % Compute new state information
    rsat_new = rsat + dr_sat;
    theta_sat_new = theta_sat + dtheta_sat;
    dr_sat_new = dr_sat + ((rsat + 0.5*dr_sat)*dtheta_sat^2 - GMs/(rsat^2)*dt^2);
    dtheta_sat_new = dtheta_sat - (2*dr_sat*dtheta_sat/(rsat + 0.5*dr_sat));
    
    % Store information for r and theta
    all_rsat = [all_rsat, rsat_new];
    all_theta_sat = [all_theta_sat, theta_sat_new];

    % Change new state information to current state information for next iteration
    rsat = rsat_new;
    dr_sat = dr_sat_new;
    theta_sat = theta_sat_new;
    dtheta_sat = dtheta_sat_new;
end;

% Plot Results
figure(1)
% Animation of all three orbits synced with time.
for k = valueoffset:200:floor(.5*length(all_theta_sat)+valueoffset)
   if k > length(all_theta_earth)
   te = all_theta_earth(abs(k-length(all_theta_earth)));
   ye = all_rearth(abs(k-length(all_theta_earth)));
   he = polar(te,ye,'-ob');
   refreshdata(he,'caller')
   drawnow;
   hold on; 
   else
   te = all_theta_earth(k);
   ye = all_rearth(k);
   he = polar(te,ye,'-ob');
   refreshdata(he,'caller')
   drawnow;
   hold on;
   end;
   
   if (k - length(all_theta_earth)) >= 0
   tv = all_theta_venus(abs(k - length(all_theta_earth)));
   yv = all_rvenus(abs(k - length(all_theta_earth)));
   hv = polar(tv,yv,'-or');
   refreshdata(hv,'caller')
   drawnow;
   hold on;
   else
   tv = all_theta_venus(length(all_theta_venus)+k-length(all_theta_earth));
   yv = all_rvenus(length(all_theta_venus)+k-length(all_theta_earth));
   hv = polar(tv,yv,'-or');
   refreshdata(hv,'caller')
   drawnow;
   hold on;
   end;
   
   ts = all_theta_sat(k-(valueoffset-1));
   ys = all_rsat(k-(valueoffset-1));
   hs = polar(ts,ys,'-ok');
   refreshdata(hs,'caller')
   drawnow;
   hold on;
  
end;

[X_earth,Y_earth] = pol2cart(all_theta_earth,all_rearth);
[X_venus,Y_venus] = pol2cart(all_theta_venus,all_rvenus);
[X_sat,Y_sat] = pol2cart(all_theta_sat,all_rsat);

figure(1)
plot(X_earth,Y_earth,'b')
hold on
plot(X_earth(40063),Y_earth(40063),'ob') %Earth Position at Sat Departure Time
hold on
plot(X_earth(8272),Y_earth(8272),'ob') %Earth Position at Sat Arrival Time
hold on
plot(X_earth(21232),Y_earth(21232),'ob') %Earth Position at Sat Arrival Time +90 days
hold on
plot(X_venus,Y_venus,'r')
hold on
plot(X_venus(19826),Y_venus(19826),'or') %Venus Position at Sat Departure Time
hold on
plot(X_venus(8272),Y_venus(8272),'or') %Venus Position at Sat Arrival Time
hold on
plot(X_venus(21232),Y_venus(21232),'or') %Venus Position at Sat Arrival Time +90 days
hold on
plot(X_sat,Y_sat,'k')
hold on
plot(X_sat(1),Y_sat(1),'ok') %Sat Position at Sat Departure Time
hold on
plot(X_sat(20801),Y_sat(20801),'ok') %Sat Position at Sat Arrival Time
axis off

figure(2)
polar(all_theta_earth,all_rearth,'b')
hold on
polar(all_theta_earth(40063),all_rearth(40063),'ob') %Earth Position at Sat Departure Time
hold on
polar(all_theta_earth(8272),all_rearth(8272),'ob') %Earth Position at Sat Arrival Time
hold on
polar(all_theta_earth(21232),all_rearth(21232),'ob') %Earth Position at Sat Arrival Time +90 days
hold on
polar(all_theta_venus,all_rvenus,'r')
hold on
polar(all_theta_venus(19826),all_rvenus(19826),'or') %Venus Position at Sat Departure Time
hold on
polar(all_theta_venus(8272),all_rvenus(8272),'or') %Venus Position at Sat Arrival Time
hold on
polar(all_theta_venus(21232),all_rvenus(21232),'or') %Venus Position at Sat Arrival Time +90 days
hold on
polar(all_theta_sat,all_rsat,'k')
hold on
polar(all_theta_sat(1),all_rsat(1),'ok') %Sat Position at Sat Departure Time
hold on
polar(all_theta_sat(20801),all_rsat(20801),'ok') %Sat Position at Sat Arrival Time



    
