% Feb 8, 2024
%
% This program runs a central force problem
%
% m x''(t) = - u(|x|) ( x / |x| )
%
% subject to x(0) and x'(0) given
%
% We take x=x(t) to be two-dimenional, and hence write an equivalent ODE
% for y=[y_1,y_2,y_3,y_4], where x'=[y_1,y_2] and x=[y_3,y_4]
%

function Euler_Newton(g,m,y_0,h,N,pause_len,skip)
%
% g,m,y_0 are for U(r)=-mg/r, so u(r)=U'(r)=mg/r^2 (Newton's Law for Planets)
% 
% pause_len and skip are parameters for plotting

  close all  % close all figure windows

  yvals = zeros(N+1,4);	% initialize a zero matrix for the y_i values
  yvals(1,:)=y_0;         % set yvals(1,:) to y_0, the initial condition

% Euler's method:
  for i=1:N
    yvals(i+1,:) = yvals(i,:) + h * f( yvals(i,:) );
  end  

% To test the accuracy, we compute the Energy at each yval(i,:)
%
  energ = zeros(N+1,1);
  for i=1:N+1
    ken_en = (1/2) * m * (yvals(i,1)^2+yvals(i,2)^2) ;
    pot_en = U( sqrt( yvals(i,3)^2 + yvals(i,4)^2 ) ) ;
    energ(i) = ken_en - pot_en ;
  end
  myplot(g,m,yvals,energ,N,h,pause_len,skip)
  Euler_Newton = yvals;
end


function ODE_funct= f1(y) 	% the function, f, in the ODE y' = f(y)
  ODE_funct = transpose( [0,0,-1,0; 0,0,0,-1; 1,0,0,0; 0,1,0,0 ] * transpose(y) );
end

function ODE_funct=f(y) 	% the function, f, in the ODE y' = f(y)
  ODE_funct = y;
  ODE_funct(3)=y(1);
  ODE_funct(4)=y(2);
  r = sqrt( y(3)^2+y(4)^2 );
  ODE_funct(1)= - u(r) * y(3) / r;
  ODE_funct(2)= - u(r) * y(4) / r;
end

function central=u(r)
  central=1/r^2;
end

function cent_int=U(r)
  cent_int=1/r;
end

function myplot(g,m,yvals,energ,N,h,pause_len,skip)
%
  for k = 1 : skip : N+1        % we "skip" over points
    pause(pause_len);
    hold off;
    fig = figure(1);
%
%           I've commented out the line below; it works better for my Mac,
%           but is machine dependent, so it's safer to use the default settings
%   fig.Position = [200 50 900 900];  % [ left bottom width height]
%
    plot(yvals(k,3),yvals(k,4),...
      'Marker','o','MarkerEdgeColor','red','MarkerFaceColor','red');
    hold on;
    plot(0,0,...
      'Marker','o','MarkerEdgeColor','black','MarkerFaceColor','black');
    hold on;
    set(gca,'xLim',[-6,6],'yLim',[-6,6]);
    tyme = sprintf('t=%5.2f, ',(k-1)*h);
    enrg0 = sprintf('E(0)=%5.2f, ',energ(1));
    enrg = sprintf('E(t)=%7.4f',energ(k));
%    title([ tyme ,enrg0 , enrg] , 'FontSize', 40);
    title([ tyme ,enrg0 , enrg] );      % this gives you the default font size
  end  
end