% 17-09-15, Numerical solution of 1-D motion where a proton is accelerated

% Initial parameters
M       = 1.7e-27;           % Mass (kg)
F       = 1e-12;             % Force (N)
c       = 3e8;               % Speed of light (m/s)

dT      = 10e-9;             % Time step (s)
EndT    = 2000e-9;           % Total time (s)

T       = 0:dT:EndT;         % Time vector (s)
Nloop   = length(T);         % Number of iterations
 
% INITIALIZE COORDINATES
% Without relativity
X    = 0;                    % (m)
Vx   = 0;                    % (m/s)
% With relativity
XR   = 0;
VxR  = 0;

% Vectors where positions and speeds are saved
XS    = zeros( Nloop, 1 );
VxS   = zeros( Nloop, 1 );

% With relativity
XRS   = zeros( Nloop, 1 );
VxRS  = zeros( Nloop, 1 );

tic

% Start itertion
CPUT  = cputime;
for I = 2:Nloop
   % Momentum for the system

   % Without relativity   
	PxNow  = M * Vx;                % Momentum: kg*m/s

   % With relativity
   Gamma  = 1 / sqrt( 1 - (VxR/c)^2 );
	PxNowR = M * VxR * Gamma;       % Momentum: kg*m/s
   
   % New momentum after time step dT
   PxNew  = PxNow  + F * dT;
   PxNewR = PxNowR + F * dT;
   
   % New speed after time step dT
   VxNew  = PxNew / M;
   VxRNew = PxNewR / M / sqrt( 1 + (PxNewR/(M*c))^2 );
   if VxRNew >= c, error('VxRNew >= c');  end

   
   % New position after time step dT
   X  = X  + 0.5* ( Vx  + VxNew  ) * dT;
   XR = XR + 0.5* ( VxR + VxRNew ) * dT;
   
   Vx   = VxNew;
   VxR  = VxRNew;

   XS(I)   = X;
   VxS(I)  = Vx;
   XRS(I)  = XR;
   VxRS(I) = VxR;
end  % for I = 2:Nloop


figure(1)
plot( T, XS,'ko-' , T,XRS,'rs--')
legend('Icke relativistiskt','Relativistiskt')
xlabel('Time (s)')
ylabel('Plats (m)')

figure(2)
plot( T, VxS,'ko-' , T,VxRS,'rs--')
legend('Icke relativistiskt','Relativistiskt')
xlabel('Time (s)')
ylabel('Hastighet (m/s)')

toc