% Probmem 1 Solution

%Here is the solution in matlab.  See the posted example from the first
%section for a handwritten solution (only with different initial
%conditions).  Hopefully this will be more helpful to you than seeing 
%another handwritten version.

clear all; close all; clc;
%initial conditions
rx0 = 1.5;
vx0 = 0.0;
ry0 = 1.0;
vy0 = 0.5;

dt = 0.1; %time interval (sec) (epsilon in problem statement)
Gm = 1;

%a) First set with xdot = b(t) (i.e. the left point on the interval)
for i = 1:3  %will calculate for three iterations
    
%Each soln will be saved in a vector i.e. rx = [rx0, rx1, rx2, rx3]
%for each iteration.  So first set the first component of the vector
%equal to the initial values
    rxa(1) = rx0;
    vxa(1) = vx0;
    rya(1) = ry0;
    vya(1) = vy0;
   %find magnitude of r at given time step
   ra(i) = sqrt(rxa(i)^2 + rya(i)^2);
   %calc vx at next time step
   vxa(i+1) = vxa(i) + dt * (-rxa(i)/(ra(i)^3));
   %calc rx at next time step
   rxa(i+1) = rxa(i) + dt * vxa(i);
   %find vy at next time step
   vya(i+1) = vya(i) + dt * (-rya(i)/(ra(i)^3));
   %calc ry at next time step
   rya(i+1) = rya(i) + dt * vya(i);
end

%b) Second set with xdot = b(t + dt) (i.e. the right point on the interval)
for i = 1:3  %will calculate for three iterations
    rxb(1) = rx0;
    vxb(1) = vx0;
    ryb(1) = ry0;
    vyb(1) = vy0;
   rb(i) = sqrt(rxb(i)^2 + ryb(i)^2);
   vxb(i+1) = vxb(i) + dt * (-rxb(i)/(rb(i)^3));
   rxb(i+1) = rxb(i) + dt * vxb(i+1);
   vyb(i+1) = vyb(i) + dt * (-ryb(i)/(rb(i)^3));
   ryb(i+1) = ryb(i) + dt * vyb(i+1);
end

%c) Third set with xdot = b((t + dt)/2) (i.e. slope at time-midpoint)
for i = 1:3  %will calculate for three iterations
    rxc(1) = rx0;
    vxc(1) = vx0;
    ryc(1) = ry0;
    vyc(1) = vy0;
   rc(i) = sqrt(rxc(i)^2 + ryc(i)^2);
   vxc(i+1) = vxc(i) + dt * (-rxc(i)/(rc(i)^3));
   %to get vx at the time midpoint, use half the time step (dt/2)
   vxc_midpt(i) = vxc(i) + dt/2 * (-rxc(i)/(rc(i)^3)); 
   rxc(i+1) = rxc(i) + dt * vxc_midpt(i);
   vyc(i+1) = vyc(i) + dt * (-ryc(i)/(rc(i)^3));
   %to get vy at the time midpoint, again use half the time step (dt/2)
   vyc_midpt(i) = vyc(i) + dt/2 * (-ryc(i)/(rc(i)^3)); 
   ryc(i+1) = ryc(i) + dt * vyc_midpt(i);
end

%d) Forth set with xdot = b((b(t)+ b(t+dt))/2) (i.e. the avg slope)
for i = 1:3  %will calculate for three iterations
    rxd(1) = rx0;
    vxd(1) = vx0;
    ryd(1) = ry0;
    vyd(1) = vy0;
   rd(i) = sqrt(rxd(i)^2 + ryd(i)^2);
   vxd(i+1) = vxd(i) + dt * (-rxd(i)/(rd(i)^3));
   vxd_avg(i) = (vxd(i) + vxd(i+1))/2; %calc avg vxc
   rxd(i+1) = rxd(i) + dt * vxd_avg(i);
   vyd(i+1) = vyd(i) + dt * (-ryd(i)/(rd(i)^3));
   vyd_avg(i) = (vyd(i) + vyd(i+1))/2; %calc avg vyc
   ryd(i+1) = ryd(i) + dt * vyd_avg(i);
end

figure()
plot(rxa,rya,'o',rxb,ryb,'o',rxc,ryc,'o',rxd,ryd,'o', 'linewidth',2)
title('Approximation of Orbital Path')
xlabel('r_x')
ylabel('r_y')
axis([1.485 1.505 0.97 1.17])
legend('xdot = b(t)','xdot = b(t+dt)','xdot = b((t+dt)/2)',...
    'xdot = b((b(t)+ b(t+dt))/2)','location', 'SouthWest')

%Display all values
rxa
rya
vxa
vya
rxb
ryb
vxb
vyb
rxc
ryc
vxc
vyc
rxd
ryd
vxd
vyd