clear all
close all
format compact


for T = 5:5:175
%for T = [15,45,90,175]


%K0 = input('alpha energy[MeV]: ');
K0 = 7.7;
Np = 200             ; % Number of points to use to plot curves
delta = 0.1          ;
z = 2                ; % Charge on the alpha
Z = 79               ; % Gold nucleus atomic number
A = 197              ; % Gold nucleus atomic mass
Malpha = 3727.3803   ; % Mass of the alpha in MeV
EM = 1.439976        ; % EM constant e^2/(4 pi epsilon 0) MeV-fm
 

%T  = input('Desired scattering angle (in degrees): ');
%T = 45;
b  = (z*Z/(2*K0))*EM*cotd(T/2);
Rc = z*Z*EM/K0; % 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*K0);

if (T > 0 & T < 180)
   kappa = (z*Z/(2*K0*b))*EM;
   rmin = b/(sqrt(1 + kappa^2) - kappa);
else
   disp('Theta must be  >0 and <180 degrees, Stopping.')
   return
end

if (rmin < (Rn + Rp))
   disp('Warning. Alpha will enter nucleus. Trajectory is approximate.')
end


disp(sprintf('\nZ target                         =  %3.0f',Z        ))
disp(sprintf(  'Z projectile                     =  %3.0f',z        ))
disp(sprintf(  'Projectile kinetic energy        =  %3.1f [MeV]',K0 ))
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(  'Impact parameter                 =  %7.2f [fm]',b   ))
disp(sprintf(  'Closest approach                 =  %7.2f [fm]',rmin))
%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 Coulomb barrier
angles = linspace(0,2*pi,Np);
xC = Rc*cos(angles);
yC = Rc*sin(angles);
plot(xC,yC,'r','linewidth',2)
s = sprintf('Nuclear surface (black), Coulomb barrier (red), K-shell radius (green)');
title(s,'fontsize',14)
xlabel('z [fm]','fontsize',16)
ylabel('\rho [fm]','fontsize',16)
hold on

% Show the nuclear surface
xN = Rn*cos(angles);
yN = Rn*sin(angles);
plot(xN,yN,'k','linewidth',2)

%Show the inner K shell orbital
xA = Ra*cos(angles);
yA = Ra*sin(angles);
plot(xA,yA,'g','linewidth',2)

uy = 800; %Our universe
ux = 1000; %Our universe
U = ux^2 + uy^2;
axis([-ux ux -uy uy])

%pause

%ZU = input('Starting z of the the projectile: ');
ZU = -ux;
disp(sprintf('Trajectory starts at (z,rho) =  (%7.2f, %7.2f)',Z,b))

Theta0 = acosd(-kappa/sqrt(kappa^2 + 1));

x = b/(ZU^2 + b^2);
theta0 = acosd((-kappa*(x + kappa) - sqrt(1 - 2*x*kappa - x^2))/(1 + kappa^2));
eps = 180 - theta0;
theta1 = T + eps;
theta = linspace(theta0,theta1,Np);
x = sind(theta) - kappa*(1 + cosd(theta));
r = b./x;
plot(r.*cosd(theta),r.*sind(theta),'k--','linewidth',2)

% This is the line from the origin to the distance of closest approach
x0 = sind(Theta0) - kappa*(1 + cosd(Theta0));
plot([0,b*cosd(Theta0)/x0],[0,b*sind(Theta0)/x0],'k--')

asymptote = 1;
if (asymptote ==1)
   % Asymptote parallel to initial and final directions, meeting up at the line
   % of closest approach
   xa = [-ux,b/tand(Theta0),b/tand(Theta0) + U*cosd(T)];
   ya = [b,b,b+U*sind(T)];
   plot(xa,ya,'k--') % Parallel with initial trajectory

end

s = sprintf( ...
  'z = %3.0f, Z = %2.0f, K = %3.1f, b = %7.3f, Theta = %7.2f',z,Z,K0,b,T);
title(s,'fontsize',16)

pause
close
clear all
end
