PhotonAttenuation(Material, Energy, Options, Thickness)

function X = PhotonAttenuation(Material, Energy, Options, Thickness)
% PhotonAttenuation provides photon attenuation in different materials.
%
% X = PhotonAttenuation(Material, Energy, Options, Thickness)
% Function providing the attenuation and energy absorption of x-ray and
% gamma-ray photons in various materials including mixtures and compounds, 
% based on NIST report 5632, by J. Hubbell and S.M. Seltzer.
%
% Input :
% 1) Material - string, number or array of strings, numbers or cells 
%     describing material type. 
%     - Element atomic number Z - in 1 to 100 range
%     - Element symbols - 'Pb', 'Fe'
%     - Element names   - 'Lead', 'Iron', 'Cesium' or 'Caesium'
%     - Some common names and full compound names - 'Water', 'Polyethylene' 
%       (see function PhysProps for more details)
%     - Compound formulas - 'H2SO4', 'C3H7NO2'- those are case sensitive 
%     - Mixtures of any of above with fractions by weight - like
%       'H(0.057444)C(0.774589)O(0.167968)' for Bakelite or  
%       'B(10)H(11)C(58)O(21)' for Borated Polyethylene (BPE-10)
% 2) Energy - Energy of the photons. Can be single energy or vector of 
%       energies. Several formats are allowed: 
%     - Energy in MeV, should be in [0.001, 20] MeV range. 
%     - Wavelengths in nano-meters. Encoded as negative numbers. The valid
%     range is 6.1982e-5 to 1.2396 nm;
%     - Continuous Spectrum - Encoded in 2 columns: column one contains 
%       energy and column two contains relative number of photons at that 
%       energy. Spectrum is assumed to be continuous and output is
%       calculated through integration using 'trapz' function.
% 3) Options - specifies what to return. String or number:
%     1 - 'mac' - function returns Mass Attenuation Coefficients in cm^2/g
%     2 - 'meac' - function returns Mass Energy-Absorption Coefficients 
%           in cm^2/g. See link below for more info:
%           http://physics.nist.gov/PhysRefData/XrayMassCoef/chap3.html
%     3 - 'cross section' or 'x' - function returns cross section in barns per
%          atom (convert to cm^2 per atom by multiplying by 10^-24). 
%          Available only for elements.
%     4 - 'mean free path' or 'mfp' - function returns mean free path (in cm) of 
%           photon in the given material . Available only for chemicals 
%           recognized by 'PhysProps' function (since density is needed).
%     5 - 'transmission' or 't' - fraction of protons absorbed by given thickness
%          of material
%     6 - 'ln_T' - log of transmission
%     7 - 'lac' - Linear Attenuation Coefficients in 1/cm same as 1/(Cross
%          Section) or mac*density.
% 4) Thickness - Thickness of material in cm. Either scalar or vector of 
%      the same length as number of materials. Negative numbers indicate
%      mass thickness measured in g/cm^2 (density*thickness). Needed only
%      if energy spectrum is used or in case of Options set to
%      'Transmission'.
% Output :
%   X  - output depends on 'Options' parameter. Columns always correspond
%        to the materials and rows correspond to energy. In case of spectrum
%        input only single row is returned.
%
% History:
%   Written by Jarek Tuszynski (SAIC), 2006
%   Inspired by John Schweppe Mathematica code available at
%     http://library.wolfram.com/infocenter/MathSource/4267/
%   Tables are based on NIST's XAAMDI and XCOM databases. See
%     PhotonAttenuationQ function for more details.
%
% Examples:
% %% Plot Cross sections of elements for different energy photons 
% figure
% Z = 1:100; % elements with Z in 1-100 range
% E = logspace(log10(0.001), log10(20), 500);  % define energy grid
% X = PhotonAttenuation(Z, E, 'cross section');
% imagesc(log10(X)); colorbar;
% title('Cross sections of elements for different energy photons');
% xlabel('Atomic Number of Elements');
% ylabel('Photon Energy in MeV');
% set(gca,'YTick',linspace(1, length(E), 10));
% set(gca,'YTickLabel',1e-3*round(1e3*logspace(log10(0.001), log10(20), 10)))
% 
% %% Plot Photon Attenuation Coefficients, using different input styles
% figure
% E = logspace(log10(0.001), log10(20), 500);  % define energy grid
% Z = {'Concrete', 'Air', 'B(10)H(11)C(58)O(21)', 100, 'Ag'};
% mac  = PhotonAttenuation(Z, E, 'mac');
% loglog(E, mac);
% legend({'Concrete', 'Air', 'BPE-10', 'Fermium', 'Silver'});
% ylabel('Attenuation in cm^2/g');
% xlabel('Photon Energy in MeV');
% title('Photon Attenuation Coefficients for different materials');

%% Process 'Energy' Input Parameter
s = size(Energy);
nEnergy = prod(s);
Spectrum = [];
if (min(s)==2)
  if (s(1)==2), Energy=Energy'; end
  Spectrum = Energy(:,2);
  Energy   = Energy(:,1);
  Spectrum = Spectrum / trapz(Energy, Spectrum); % normalize spectrum
  nEnergy  = 1; % this is a spectrum output will be an array instead of matrix
end
Energy  = Energy(:);
if (any(Energy<0)) 
  idx = find(Energy<0);
  PlanckConstant = 4.1350e-021; % Planks constant in MeV*sec
  SpeedOfLight   = 299792458E9; % in nano meters per second
  WaveLength     = -Energy(idx);
  Energy(idx)    = (PlanckConstant*SpeedOfLight) ./ WaveLength;
end
if (any(Energy<0.0009) || any(Energy>21))
  warning('PhotonAttenuation:wrongEnergy',...
    'Warning in PhotonAttenuation function: energy is outside of the recomended range from 0.001 MeV to 20 MeV');
end

%% Process 'Material' Input Parameter
if (ischar   (Material)), Material = {Material}; end
if (isnumeric(Material)), Material = num2cell(Material); end
nName = length(Material);

%% Process 'Options' Input Parameter
if (nargin<3), Options = 1; end
if (ischar(Options))
  X = {'mac','meac', 'cross section', 'mean free path', 'transmission', 'ln_t', 'lac'...
       'mac','meac', 'x', 'mfp', 't', 'log_t', 'lac'};
  Options = find(strcmpi( Options, X),1);
  if (Options>7), Options=Options-7; end
else
  if (Options>7), Options=[]; end
end
if (isempty(Options)) 
  error(['Error in PhotonAttenuation function: Options parameter was not recognized: ', Options]);
end
saveMac = (Options==2)+1;

%% Process 'Thickness' Input Parameter
if (nargin<4), Thickness=1; end
if (length(Thickness)==1 && nName>1)
  Thickness = repmat(Thickness(1), 1, nName);
end
if (length(Thickness)>1 && nName==1)
  nName = length(Thickness);
  Material = repmat(Material(1), 1, nName);
end
if (length(Thickness)~=nName)
  error('Error in PhotonAttenuation function: missmatch between lengths or Material and Thickness arrays');
end


%% Initialize variables
if (isempty(Spectrum)), X = zeros(nEnergy,nName);
else                    X = zeros(1      ,nName); end
u = 1.6605402E-24 ; % atomic mass unit (1/12 of the mass of C-12)
barns = 1E24;      % barns 

%% Main Loop
Atoms = [];
for i = 1:nName
  % Parse name
  if (~isempty(Atoms) && length(Atoms)==length(Material{i}) && all(Atoms==Material{i}))
    if (Options<3 && isempty(Spectrum)) % a shortcut for common configuration
      X(:,i) = mu; 
      continue;
    end
  end
  
  Atoms   = Material{i};
  Z       = Atoms;
  Ratios  = 1;
  if (isempty(Z))
    error('Error in PhotonAttenuation function: compound formula was empty');
  end
  if (ischar(Atoms)) % check if this is element name or known compound
    [Z Ratios] = ParseChemicalFormula(Atoms);
    if (isempty(Z))
      error(['Error in PhotonAttenuation function: compound formula was not recognized: ', Atoms]);
    end
    MatStr = Atoms;
  else
    MatStr = num2str(Atoms);
  end
  if (length(Z)>1 && Options==3)
    error('Error in PhotonAttenuation function: Cross section "Option" is only supported for elements.');
  end

  % Look up the data, average contributions of different elements and save
  mu = PhotonAttenuationQ(Z, Energy, saveMac);  
  Ratios = Ratios(:) / sum(Ratios); % normalize ratios so they add up to one
  % if mu is 2D (spectrum and compound) it will become 1D 
  % if mu is 1D (because of compound) it will become scalar
  % if mu is 1D (because of spectrum) it will not change
  % if mu is a scalar (monoenergetic & element) it will not change
  mu  = mu * Ratios; 
  if (Options<3 && isempty(Spectrum)) % a shortcut for common configuration
    X(:,i) = mu; 
    continue;
  end
  
  % if desired output needs density than try to look it up
  PP      = PhysProps(Atoms);
  Density = PP{1,2}; 
  A       = Z/PP{1,1}; % atomic molar mass in g/mol
 
  % if Energy has a spectrum than integrate
  if(Thickness(i)>0), MassThick =  Thickness(i)*Density;
  else                MassThick = -Thickness(i); end % user provided linear thickness
  MassThick = max(MassThick, 1e-6);
  T = exp(-mu*MassThick); % always convert to transmission
  mac = mu;
  if (~isempty(Spectrum))
    T  = trapz(Energy, T.*Spectrum); % integrate over energy (T become a scalar)
    mac = -log(T)/MassThick;
  end
  
  % cast in chosen format
  switch (Options)
    case {1,2} % mac & meac
      X(:,i) = mac;
    case 3 % 'cross section'
      X(:,i) = mac * (barns*u*A);
    case 4 % 'mean free path'
      X(:,i) = 1./(mac*Density);
    case 5 % 'transmission'
      X(:,i) = T;
    case 6 % 'log T'
      X(:,i) = -log(T);
    case 7 % 'linear attenuation coefficiant'
      X(:,i) = mac*Density;
  end
  
  if (any(isnan(X(:,i))))
    error(['Error in PhotonAttenuation function: Physical properties of "',...
           MatStr,'" not found in PhysProps function.']);
  end

end