% randomWalks.m A Matlab script M-file demonstrating many type of random walks
% 
% FileName      : randomWalks.m
% Creation      : 09/29/2016
% Modifications : September 29, 2016 Initial creation
% Author(s)     : Alex Bielajew
%

clear all
close all
format long
format compact
clc

disp(sprintf('Several types of random walks are performed'))
disp(sprintf('==========================================='))
disp(sprintf('1: Equal steps, 2D, isotropic'))
disp(sprintf('2: Equal steps, 3D, semi-isotropic'))            

Sim = input('Which simulation [1:4] would you like to see? ');

if ((floor(Sim) ~= Sim) || (Sim < 1) || (Sim > 4))
   disp('Input must be integer an between 0 and 4 incusive. Try again.')
   return
end

clc

if (Sim == 1)
   disp(sprintf('1: Equal steps, 2D, isotropic'))
   disp(sprintf('============================='))
   Nh   = input('Number of histories: ');
   smax = input('Number of step/history: ');

   % Size the trajectory arrays
   x = zeros(1,smax+1); % smax steps => smax+1 endpoints
   y = x;
   
   if (Nh <= smax)
      figure(1)
      hold on
   end
   
   %Loop over the number of histories
   maxmax = 0; % Used for axis scaling

   tic
   for n = 1:Nh
   
      if (~mod(n,10000))  % Throbber to thow that it has not crashed
         t = toc;
         tic;
         disp(sprintf('After %g histories, time = %g s',n,t))
      end

      x(1) = 0;
      y(1) = 0;
      
      % Loop over the individual steps/history
      for i = 1:smax
         [u,v] = azimuthal; % "azimuthal" is in the code library
         % Note that the first step is deflected before transport
         x(i+1) = x(i) + u;
         y(i+1) = y(i) + v;
      end
      
      % Display the histories cumulatively
      if (Nh <= smax)
         plot(x,y,'k-o')
         maxx = max(abs(x));
         maxy = max(abs(y));
         maxmax = max(maxmax,max(maxx, maxy));
         axis equal
         axis([-maxmax, maxmax, -maxmax, maxmax])
         if (n == 1)
            title(sprintf(' After %g history',n))
         else
            title(sprintf(' After %g histories',n))
         end
         if (n <= 10)
            pause
         else
            pause(0.01)
         end
         %close % Uncomment this to see individual histories
      end
      
      xend(n) = x(end); % Used for ploting the endpoints only
      yend(n) = y(end);
      
   end % End of loop over histories
   
   choice = input(' Plot all the endpoints? (0 => no) ');
   if (choice)
      figure(2)
      plot(xend,yend,'ko')
   end
   
   % Calculate some stats
   xbar   = sum(xend)/Nh;
   x2bar  = sum(xend.^2)/Nh;
   s2x    = (x2bar - xbar^2);
   s2xbar = (x2bar - xbar^2)/Nh;
   x3bar  = sum(xend.^3)/Nh; 
   x4bar  = sum(xend.^4)/Nh;
   ybar   = sum(yend)/Nh;
   y2bar  = sum(yend.^2)/Nh;
   y3bar  = sum(yend.^3)/Nh ;
   y4bar  = sum(yend.^4)/Nh;
   s2y    = (y2bar - ybar^2);
   s2ybar = (y2bar - ybar^2)/Nh;
   disp(' ')
   disp(sprintf('<x>  = %g +/- %g',xbar, sqrt(s2xbar)))
   disp(sprintf('<y>  = %g +/- %g',ybar, sqrt(s2ybar)))
   disp(' ')
   disp(sprintf('<x^2>  = %g',x2bar))
   disp(sprintf('<y^2>  = %g',y2bar))
   disp(' ')
   disp(sprintf('<x^3>  = %g',x3bar))
   disp(sprintf('<y^3>  = %g',y3bar))
   disp(' ')
   disp(sprintf('<x^4>  = %g',x4bar))
   disp(sprintf('<y^4>  = %g',y4bar))
   disp(' ')

elseif (Sim == 2)

   disp(sprintf('1: Equal steps, 3D, semi.isotropic'))
   disp(sprintf('=================================='))
   Nh   = input('Number of histories: ');
   smax = input('Number of step/history: ');

   % Size the trajectory arrays
   x = zeros(1,smax+1); % smax steps => smax+1 endpoints
   y = x;
   z = y;
   
   if (Nh <= smax)
      figure(1)
      hold on
   end
   
   constraint = input('Maximum scattering angle (in degrees): ');
   maxTheta = constraint*pi/180; % Maximum angle in radians
   cosmaxTheta = cos(maxTheta);

   %Loop over the number of histories

   tic
   for n = 1:Nh
   
      if (~mod(n,10000))  % Throbber to thow that it has not crashed
         t = toc;
         tic;
         disp(sprintf('After %g histories, time = %g s',n,t))
      end

      % Take the first step
      x(1) = 0;
      y(1) = 0;
      z(1) = 0;
      x(2) = 0;
      y(2) = 0;
      z(2) = 1;
      u = [0,0,1]; % Initial direction vector
      
      % Loop over the individual steps/history
      for i = 2:smax
         ws = 1 - rand*(1 - cosmaxTheta); % This is semi-isotropic scattering
         u = rotate(u,ws); % "rotate" is in the code library
         unit = sum(u.*u);
         if (abs(1-unit) >  10^(-14)) % Norm check
            unit - 1
         end
         % In this scheme, deflection is after the displacement
         x(i+1) = x(i) + u(1);
         y(i+1) = y(i) + u(2);
         z(i+1) = z(i) + u(3);
      end
      
      % Display the histories cumulatively
      if (Nh <= smax)
         plot3(z,x,y,'k-o')
         if (n == 1)
            title(sprintf('After %g history',n))
         else
            title(sprintf('After %g histories',n))
         end
         xlabel('z')
         ylabel('x')
         zlabel('y')
         if (n <= 10)
            pause
         else
            pause(0.01)
         end
         %close % Uncomment this to see individual histories
      end
      
      xend(n) = x(end); % Used for ploting the endpoints only
      yend(n) = y(end);
      zend(n) = z(end);
      
   end % End of loop over histories
   
   choice = input('Plot all the endpoints? (0 => no) ');
   if (choice)
      figure(2)
      plot3(zend,xend,yend,'ko')
      xlabel('z')
      ylabel('x')
      zlabel('y')
   end
   
   % Calculate some stats
   xbar = sum(xend)/Nh/smax;
   x2bar = sum(xend.^2)/Nh/smax;
   x3bar = sum(xend.^3)/Nh/smax; 
   x4bar = sum(xend.^4)/Nh/smax/(smax - 1);
   ybar  = sum(yend)/Nh/smax;
   y2bar = sum(yend.^2)/Nh/smax;
   y3bar = sum(yend.^3)/Nh/smax ;
   y4bar = sum(yend.^4)/Nh/smax/(smax - 1);
   zbar  = sum(zend)/Nh/smax;
   z2bar = sum(zend.^2)/Nh/smax;
   z3bar = sum(zend.^3)/Nh/smax;
   z4bar = sum(zend.^4)/Nh/smax/(smax - 1);
   disp(' ')
   disp(sprintf('<x>  = %g',xbar))
   disp(sprintf('<y>  = %g',ybar))
   disp(sprintf('<z>  = %g',zbar))
   disp(' ')
   disp(sprintf('<x^2>  = %g',x2bar))
   disp(sprintf('<y^2>  = %g',y2bar))
   disp(sprintf('<z^2>  = %g',z2bar))
   disp(' ')
   disp(sprintf('<x^3>  = %g',x3bar))
   disp(sprintf('<y^3>  = %g',y3bar))
   disp(sprintf('<z^3>  = %g',z3bar))
   disp(' ')
   disp(sprintf('<x^4>  = %g',x4bar))
   disp(sprintf('<y^4>  = %g',y4bar))
   disp(sprintf('<z^4>  = %g',z4bar))
   disp(' ')
  
elseif (Sim == 3)

elseif (Sim == 4)

   
end