%
%   PTSRC2G.M  Plots the flux profiles associated with a Point Source of  
%              Fast Neutrons in an Infinite Moderator
%
%   This file plots the fast and thermal flux profiles associated with a point 
%   source of fast neutrons diffusing in an infinite moderating medium. Data 
%   for water, graphite, and beryllium are available (from Lamarsh) for comparison 
%   of how the diffusion parameters of different materials affect the diffusion and
%   attenuation of neutrons.
%
%   A plot of the fast-to-thermal flux ratio is also generated.
%
%   File prepared by J. R. White, UMass-Lowell (Jan. 2015)
%

      clear all,  close all,  nfig = 0;
% 
%   define material properties with data from Lamarsh 3rd Ed.(Tables 5.2 and 5.3)
%   --> order: water, graphite, and beryllium
      DD1 = [1.13 1.016 0.562];  % fast diffusion coeff (cm)
      TT = [27 368 102];         % thermal neutron age (cm^2)
      TTSR = sqrt(TT);           % square root of neutron age (cm)
      DD2 = [0.16 0.84 0.50];    % thermal diffusion coeff (cm)
      LLS = [8.1 3500 480];      % diffusion area (cm^2)
      LL = sqrt(LLS);            % diffusion length (cm)
      TITL = ['Water    '; 'Graphite '; 'Beryllium'];
%
%   select moderator material
      K = menu('Choose Moderator Matl', ...
               '           Water           ', ...
               '         Graphite          ', ...
               '         Beryllium         ');
      D1 = DD1(K);  T = TT(K);    TSR = TTSR(K);
      D2 = DD2(K);  LS = LLS(K);  L = LL(K);   titl = TITL(K,:);
%
%   select type of plot
      IPLT = menu('Type of Flux Plot?', '   Linear Y-Axis   ','   Log Y-Axis   ');
%
%   define radial grid  (use R approx 5*fast diffusion length)
%   Note: start at about 1 cm to avoid the singularity at r = 0   
      N = 200;   R = ceil(5*TSR/10)*10;    r = linspace(1,R,N);
%
%   compute fluxes and fast to thermal flux ratio (for unit source strength)
      c1 = 1/(4*pi*D1);   c2 = LS/(4*pi*D2*(LS-T));
      phi1 = c1*exp(-r/TSR)./r;
      phi2 = c2*(exp(-r/L) - exp(-r/TSR))./r;
      ratio = phi1./phi2;
%
%   plot flux profiles
      nfig = nfig+1;  figure(nfig)
      if  IPLT == 1, 
        subplot(2,1,1),plot(r,phi1,'r-',r,phi2,'g--','LineWidth',2),grid
      end
      if  IPLT == 2, 
        subplot(2,1,1),semilogy(r,phi1,'r-',r,phi2,'g--','LineWidth',2),grid
      end
      title(['Flux & Ratio Profiles for a Point Source of Fast Neutrons in ',titl])
      ylabel('Flux per Unit Source '),legend('fast','thermal')
%
%   plot fast-to-thermal flux ratio
      subplot(2,1,2),plot(r,ratio,'b-','LineWidth',2),grid
      xlabel('Radial Distance from Source (cm)'),ylabel('\phi_1/\phi_T')
%
%   end of problem