clear all
close all
format compact

Nt = 10              ; % Number of trajectories to plot
Np = 200             ; % Number of points to use to plot curves
Rbeam = 1            ; % Radius of the beam
delta = 0.1          ;
Z = 79               ; % Gold nucleus atomic number
A = 197              ; % Gold nucleus atomic mass
z = 2                ; % alpha-particle atomic number
Malpha = 3727.3803   ; % Mass of the alpha in MeV
EM = 1.439976        ; % EM constant e^2/(4 pi epsilon0)
Rdet = 100           ; % Radius of the detector
 

K  = input('alpha energy[MeV]: ');
Rc = z*Z*EM/K; % Coulomb barrier, used as dimension scale unit
Rn = 1.2*(A^(1/3)); % Nuclear radius
Rp = 1.2*(4^(1/3)); % Projectile radius
Ra = 52900/Z; % K electron radius, 52900 is the Bohr radius in fm
dB = 1239.853/sqrt(2*3727.409*K);

disp(sprintf('\nZ target                         =  %3.0f',Z))
disp(sprintf(  'Z projectile                     =  %3.0f',z))
disp(sprintf(  'Projectile kinetic energy        =  %3.0f [MeV]',K))
disp(sprintf(  'Projectile de Broglie wavelength =  %7.2f [fm]',dB))
disp(sprintf(  'Projectile radius                =  %7.2f [fm]',Rp))
disp(sprintf(  'Coulomb barrier                  =  %7.2f [fm]',Rc))
disp(sprintf(  'Nuclear radius                   =  %7.2f [fm]',Rn))
disp(sprintf(  'K-shell radius                   =  %7.2f [fm]',Ra))
disp(sprintf('Note. From now on all distances to be measured in units of the Coulomb barrier.\n\n'))

Equal = z*Z*EM/(Rn);

% Show the cross section
angles = linspace(pi/Np,pi,Np); % Avoids zero angle
dS = ((Rc/4)^2)./(sin(angles/2).^4);
plot(180*angles/pi,log10(dS))
s = sprintf('Rutherford differential cross section');
title(s,'fontsize',16)
s = sprintf('Scattering angle [degrees]');
xlabel(s,'fontsize',16)
s = sprintf('Common log of cross section [fm^2]');
ylabel(s,'fontsize',16)
pause

% Show the Coulomb barrier
angles = linspace(0,2*pi,Np);
xC = cos(angles);
yC = sin(angles);
plot(xC,yC,'r','linewidth',2)
s = sprintf('Scaling unit, 1 =  %7.3f fm',Rc);
title(s,'fontsize',16)
s = sprintf('Nuclear surface (black), Coulomb barrier (red), K-shell radius (green)');
xlabel(s,'fontsize',16)
hold on

% Show the nuclear surface
scale = Rn/Rc;
plot(scale*xC,scale*yC,'k','linewidth',2)

%Show the inner K shell orbital
scale = Ra/Rc; % 52900 is the Bohr radius in fm
plot(scale*xC,scale*yC,'g','linewidth',2)
%axis([-scale,scale,-scale,scale])
pause 

R = input('Longitudinal distance from target to start trajectory: ');
disp(sprintf('Trajectory starts at longitudinal distance =  %7.2f',R))

maxD = 0;
minD = -3;
ND = 100;

for r0 = R*logspace(maxD,minD,ND)
   delta = atan(r0/R);
   Theta0 = pi - delta; % Starting angle
   ct1 = 1 + cos(Theta0);
   st = sin(Theta0);
   rad = st^2 - 2*ct1/sqrt(R^2 + r0^2);
   if (rad < 0)
      disp(' Perpendicular distance from the z-axis too small, try again.')
      b = 1;
   else
      b = ct1/(st - sqrt(rad));
      ThetaScat = 2*acot(2*b);
      t = linspace(Theta0, ThetaScat + delta/2,Np);
      r = 1./(sin(t)/b - (1 + cos(t))/2/b^2);
      hold on
      plot(r.*cos(t),r.*sin(t))
      axis([-R,R,0,R+2])
      s = sprintf('z (scaled)');
      xlabel(s,'fontsize',16)
      s = sprintf('cylindrical distance (scaled)');
      ylabel(s,'fontsize',16)
      pause(0.5)
      s = sprintf( ...
         'z = %3.0f, Z = %2.0f, K = %3.0f, b = %7.3f, Theta = %7.2f',z,Z,K,b,ThetaScat*180/pi);
      title(s,'fontsize',16)
   end
end

% Plot the envelope
b = R*logspace(maxD,minD,ND);
t = 2*atan(1./b);
r = 1./(sin(t)./b - (1 + cos(t))./b./b/2);
plot(r.*cos(t), r.*sin(t),'g','linewidth',2)

pause
close(1)

% 3D plot

meshsize = 30;

% Show the Coulomb barrier
[x,y,z] = sphere(meshsize);
mesh(x,y,z);
hold on


maxD = 0;
minD = -3;
ND = 100;

while(1)
   r0 = R*rand + R*10^minD;
   delta = atan(r0/R);
   Theta0 = pi - delta; % Starting angle
   ct1 = 1 + cos(Theta0);
   st = sin(Theta0);
   rad = st^2 - 2*ct1/sqrt(R^2 + r0^2);
   if (rad < 0)
      disp(' Perpendicular distance from the z-axis too small, try again.')
   else
      phi = 2*rand*pi;
      b = ct1/(st - sqrt(rad));
      ThetaScat = 2*acot(2*b);
      if (ThetaScat >= pi/2)
         linecolor = 'r';
      elseif (ThetaScat >= pi/4)
         linecolor = 'g';
      elseif (ThetaScat >= pi/36)
         linecolor = 'b';
      else
         linecolor = 'k';
      end
      t = linspace(Theta0, ThetaScat + delta/2,Np);
      r = 1./(sin(t)/b - (1 + cos(t))/2/b^2);
      hold on
      plot3(r.*sin(t)*cos(phi), r.*sin(t)*sin(phi), r.*cos(t),linecolor)
      axis([-R-2,R+2,-R-2,R+2,-R-2,R+2])
      s = sprintf('x (scaled)');
      xlabel(s,'fontsize',16)
      s = sprintf('y (scaled)');
      ylabel(s,'fontsize',16)
      s = sprintf('z (scaled)');
      zlabel(s,'fontsize',16)
      pause(0.01)
      s = sprintf( ...
         'z = %3.0f, Z = %2.0f, K = %3.0f, b = %7.3f, Theta = %7.2f',z,Z,K,b,ThetaScat*180/pi);
      title(s,'fontsize',16)
   end
end
