% recursion for the Kepler problem with drag force (page32)
k=1
m=1  
dt=0.02
t_max=50
F_x=inline('-k*x/(x^2+y^2)^1.5','x','y','k');
F_y=inline('-k*y/(x^2+y^2)^1.5','x','y','k');
% generalize the drag force to allow dependence on powers of v
% this requires two separate components
Fdr_x=inline('-0.01*vx','vx','vy');
Fdr_y=inline('-0.01*vy','vx','vy');
%Fdr_x=inline('-0.01*vx*sqrt(vx^2+vy^2)^3','vx','vy');
%Fdr_y=inline('-0.01*vy*sqrt(vx^2+vy^2)^3','vx','vy');
N=floor(t_max/dt)
xV(1)=1
xV(2)=1
yV(1)=0
yV(2)=1.15*dt
for j=2 : N-1
    vx=(xV(j)-xV(j-1))/dt;
    vy=(yV(j)-yV(j-1))/dt;
    xV(j+1)=2*xV(j)-xV(j-1)+dt^2/m*( F_x(xV(j),yV(j),k) + Fdr_x(vx,vy) );
    yV(j+1)=2*yV(j)-yV(j-1)+dt^2/m*( F_y(xV(j),yV(j),k) + Fdr_y(vx,vy) );
end
for j=1 : N
    if floor(j/10)*10==j
        xP(j/10)=xV(j);
        yP(j/10)=yV(j);
    end
end
    
plot(xP,yP,'b+'),title('Stroboscopic trajectory plot'),figure
% the figure switch directs Matlab to keep this figure window
% Note that the plot windows will be stacked on top of each other.

% now define expressions for kinetic and potential energies:
for j=2 : N-1
  t(j)=(j-1)*dt;
  r=sqrt(xV(j)^2+yV(j)^2);  
  v_x=(xV(j+1)-xV(j-1))/(2*dt);
  v_y=(yV(j+1)-yV(j-1))/(2*dt);
  T=0.5*m*(v_x^2+v_y^2);
  V=-k/r;
  E_m(j)=T+V;
  KE(j)=T;
  L_z(j)=m*(xV(j)*v_y-yV(j)*v_x);
  Ecc(j)=sqrt(1+2*E_m(j)*L_z(j)^2/m/k^2);
end
plot(t,[E_m; KE],'.'),
title('Kinetic Energy (top), Total Energy (bottom) versus Time'),figure

plot(t,L_z,'r.'),title('Angular Momentum versus Time'),figure
% the next plot shows whether the trajectory loses ellipticity:
% a circular orbit has zero eccentricity.
plot(t,Ecc,'r.'),title('Eccentricity versus Time')