diff --git a/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8d.xml b/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8d.xml
new file mode 100644
index 0000000..7a6326b
--- /dev/null
+++ b/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8d.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
\ No newline at end of file
diff --git a/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8p.xml b/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8p.xml
new file mode 100644
index 0000000..b188e5a
--- /dev/null
+++ b/resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8p.xml
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/src/materials/cross_section.m b/src/materials/cross_section.m
new file mode 100644
index 0000000..4d002b3
--- /dev/null
+++ b/src/materials/cross_section.m
@@ -0,0 +1,37 @@
+function cs = cross_section(Z, E) % Once tested change this function to evaluate the cross section for a vector of energies
+ % CROSS_SECTION Get the cross section of the material for a given energy (CREDIT: Geant4)
+ % The values, formulae and code is taken directly from
+ % https://github.com/Geant4/geant4/blob/master/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
+ if E < 0.1; cs = 0; return; end % Below 100 eV, we are beyond the limit of the cross section table -> 0
+ % See https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsReferenceManual/html/electromagnetic/gamma_incident/compton/compton.html
+
+ persistent a b c d1 d2 d3 d4 e1 e2 e3 e4 f1 f2 f3 f4
+ if isempty(a) % Initialize the constants if they are not already
+ a = 20.0; b = 230.0; c = 440.0; % Unitless
+ d1= 2.7965e-25; d2=-1.8300e-25; % cm^2 (e-24 for the barn)
+ d3= 6.7527e-24; d4=-1.9798e-23; % cm^2 (e-24 for the barn)
+ e1= 1.9756e-29; e2=-1.0205e-26; % cm^2 (e-24 for the barn)
+ e3=-7.3913e-26; e4= 2.7079e-26; % cm^2 (e-24 for the barn)
+ f1=-3.9178e-31; f2= 6.8241e-29; % cm^2 (e-24 for the barn)
+ f3= 6.0480e-29; f4= 3.0274e-28; % cm^2 (e-24 for the barn)
+ end
+ T0 = zeros(size(Z)) + 40; % Special case for hydrogen (KeV)
+ T0(Z > 1.5) = 15; % KeV
+
+ X = max(E, T0) ./ constants.em_ee; % Unitless
+ p1Z = Z.*(d1 + e1.*Z + f1.*Z.*Z); p2Z = Z.*(d2 + e2.*Z + f2.*Z.*Z); % cm^2
+ p3Z = Z.*(d3 + e3.*Z + f3.*Z.*Z); p4Z = Z.*(d4 + e4.*Z + f4.*Z.*Z); % cm^2
+
+ cs = p1Z.*log(1.+2.*X)./X + (p2Z + p3Z.*X + p4Z.*X.*X)./(1. + a.*X + b.*X.*X + c.*X.*X.*X); % cm^2
+
+ if E < T0
+ X = (T0+1) ./ constants.em_ee; % Unitless
+ sigma = p1Z.*log(1.+2.*X)./X + (p2Z + p3Z.*X + p4Z.*X.*X)./(1. + a.*X + b.*X.*X + c.*X.*X.*X); % cm^2
+ c1 = -T0.*(sigma-cs)./cs; % Unitless
+ c2 = zeros(size(Z)) + 0.150;
+ c2(Z > 1.5) = 0.375-0.0556.*log(Z(Z>1.5));
+
+ y = log(E./T0); % Unitless
+ cs = cs .* exp(-y.*(c1+c2.*y)); % cm^2
+ end
+end
\ No newline at end of file
diff --git a/src/materials/material_attenuation.m b/src/materials/material_attenuation.m
index 9fe3551..ecc7148 100644
--- a/src/materials/material_attenuation.m
+++ b/src/materials/material_attenuation.m
@@ -7,7 +7,10 @@
mass_fractions % Mass fractions of the elements in the material
density % Density of the material in g/cm^3
mu_from_energy % A function handle to get the linear attenuation coefficient from energy. Only used if use_mex is false
- use_mex % A flag to indicate if the MEX implementation of the photon_attenuation package is available
+ end
+
+ properties (Access=private, Constant)
+ use_pa_mex = ~~exist('photon_attenuation_mex', 'file');% A flag to indicate if the MEX implementation of the photon_attenuation package is available
end
properties (Access=private, Constant)
@@ -76,15 +79,14 @@
else
error('MATLAB:invalidInput', 'Invalid number of input arguments. Use either material("material_name") or material("material_name", atomic_numbers, mass_fractions, density).');
end
- self.use_mex = ~~exist('photon_attenuation_mex', 'file'); % Use the MEX implementation of the photon_attenuation package if available
- if ~self.use_mex
+ if ~self.use_pa_mex
self.mu_from_energy = get_photon_attenuation(self.atomic_numbers, self.mass_fractions, self.density);
end
end
function mu = get_mu(self, energy)
% GET_MU Get the linear attenuation coefficient of the material for a given energy
- if self.use_mex
+ if self.use_pa_mex
mu = photon_attenuation_mex(self.atomic_numbers, self.mass_fractions, self.density, energy);
else
mu = self.mu_from_energy(energy);
@@ -93,49 +95,9 @@
function mfp = mean_free_path(self, E)
% MEAN_FREE_PATH Get the mean free path of the material for a given energy
- mfp = 1 / (constants.N_A * self.density * sum(self.mass_fractions .* ...
- material_attenuation.cross_section(self.atomic_numbers, E) ./ ...
- self.atomic_masses(self.atomic_numbers)));
- end
- end
-
- methods (Static)
- function cs = cross_section(Z, E) % Once tested change this function to evaluate the cross section for a vector of energies
- % CROSS_SECTION Get the cross section of the material for a given energy (CREDIT: Geant4)
- % The values, formulae and code is taken directly from
- % https://github.com/Geant4/geant4/blob/master/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
- if E < 0.1; cs = 0; return; end % Below 100 eV, we are beyond the limit of the cross section table -> 0
- % See https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsReferenceManual/html/electromagnetic/gamma_incident/compton/compton.html
-
- persistent a b c d1 d2 d3 d4 e1 e2 e3 e4 f1 f2 f3 f4
- if isempty(a) % Initialize the constants if they are not already
- a = 20.0; b = 230.0; c = 440.0; % Unitless
- d1= 2.7965e-25; d2=-1.8300e-25; % cm^2 (e-24 for the barn)
- d3= 6.7527e-24; d4=-1.9798e-23; % cm^2 (e-24 for the barn)
- e1= 1.9756e-29; e2=-1.0205e-26; % cm^2 (e-24 for the barn)
- e3=-7.3913e-26; e4= 2.7079e-26; % cm^2 (e-24 for the barn)
- f1=-3.9178e-31; f2= 6.8241e-29; % cm^2 (e-24 for the barn)
- f3= 6.0480e-29; f4= 3.0274e-28; % cm^2 (e-24 for the barn)
- end
- T0(Z < 1.5) = 40; % Special case for hydrogen (KeV)
- T0(Z > 1.5) = 15; % KeV
-
- X = max(E, T0) ./ constants.em_ee; % Unitless
- p1Z = Z.*(d1 + e1.*Z + f1.*Z.*Z); p2Z = Z.*(d2 + e2.*Z + f2.*Z.*Z); % cm^2
- p3Z = Z.*(d3 + e3.*Z + f3.*Z.*Z); p4Z = Z.*(d4 + e4.*Z + f4.*Z.*Z); % cm^2
-
- cs = p1Z.*log(1.+2.*X)./X + (p2Z + p3Z.*X + p4Z.*X.*X)./(1. + a.*X + b.*X.*X + c.*X.*X.*X); % cm^2
-
- if E < T0
- X = (T0+1) ./ constants.em_ee; % Unitless
- sigma = p1Z.*log(1.+2.*X)./X + (p2Z + p3Z.*X + p4Z.*X.*X)./(1. + a.*X + b.*X.*X + c.*X.*X.*X); % cm^2
- c1 = -T0.*(sigma-cs)./cs; % Unitless
- if Z > 1.5; c2 = 0.375-0.0556.*log(Z); %Unitless
- else; c2 = 0.150; %Unitless
- end
- y = log(E./T0); % Unitless
- cs = cs .* exp(-y.*(c1+c2.*y)); % cm^2
- end
+ mfp = 1 / (constants.N_A * self.density * ...
+ sum(self.mass_fractions .* cross_section(self.atomic_numbers, E) ...
+ ./ self.atomic_masses(self.atomic_numbers)));
end
end
end
\ No newline at end of file
diff --git a/tests/scatter_tests.m b/tests/scatter_tests.m
index c2b18ab..7d6bece 100644
--- a/tests/scatter_tests.m
+++ b/tests/scatter_tests.m
@@ -30,30 +30,29 @@ function test_cross_section(tc)
tol = 0.11; % https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsReferenceManual/html/electromagnetic/gamma_incident/compton/compton.html#id276
else
tol = 0.06; % https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsReferenceManual/html/electromagnetic/gamma_incident/compton/compton.html#id276
- htol = tol;
end
- cs = material_attenuation.cross_section(Z(2), E(i));
+ cs = cross_section(Z(2), E(i));
tc.verifyEqual(cs, carbon_cs(i), 'RelTol', tol);
mfp = carbon.mean_free_path(E(i));
c_exp = 12.011 / (carbon_cs(i)*constants.N_A);
tc.verifyEqual(mfp, c_exp, 'RelTol', tol);
- cs = material_attenuation.cross_section(Z(3), E(i));
+ cs = cross_section(Z(3), E(i));
tc.verifyEqual(cs, silicon_cs(i), 'RelTol', tol);
mfp = silicon.mean_free_path(E(i));
s_exp = 28.085 / (silicon_cs(i)*constants.N_A*0.4);
tc.verifyEqual(mfp, s_exp, 'RelTol', tol);
- cs = material_attenuation.cross_section(Z(4), E(i));
+ cs = cross_section(Z(4), E(i));
tc.verifyEqual(cs, iron_cs(i), 'RelTol', tol);
mfp = iron.mean_free_path(E(i));
i_exp = 55.845 / (iron_cs(i)*constants.N_A*0.3);
tc.verifyEqual(mfp, i_exp, 'RelTol', tol);
- cs = material_attenuation.cross_section(Z(5), E(i));
+ cs = cross_section(Z(5), E(i));
tc.verifyEqual(cs, iodine_cs(i), 'RelTol', tol);
mfp = iodine.mean_free_path(E(i));
@@ -61,12 +60,12 @@ function test_cross_section(tc)
tc.verifyEqual(mfp, io_exp, 'RelTol', tol);
if i > 1 % This is the far end of the fit (lowest z and energy) so it is the worst - so bad testing would have pointlessly high errors
- cs = material_attenuation.cross_section(Z(1), E(i));
- tc.verifyEqual(cs, hydrogen_cs(i), 'RelTol', htol);
+ cs = cross_section(Z(1), E(i));
+ tc.verifyEqual(cs, hydrogen_cs(i), 'RelTol', tol);
mfp = hydrogen.mean_free_path(E(i));
h_exp = 1.008 / (hydrogen_cs(i)*constants.N_A);
- tc.verifyEqual(mfp, h_exp, 'RelTol', htol);
+ tc.verifyEqual(mfp, h_exp, 'RelTol', tol);
cs = comined.mean_free_path(E(i));
comb_exp = 1/(0.7*(0.1/(h_exp) + 0.3/(c_exp) + 0.2/(s_exp*0.4) + 0.3/(i_exp*0.3) + 0.1/(io_exp)));