Separate cross section function from class

resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8d.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design"/>
+    </Category>
+</Info>

resources/project/FjfLEenEuSm4QUSLnwd8ncPq2CI/U_LqJymK3W1-yahSVyqiNLVIya8p.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info location="cross_section.m" type="File"/>

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

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

tests/scatter_tests.m

@@ -30,43 +30,42 @@ 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));
                 io_exp = 126.90447 / (iodine_cs(i)*constants.N_A);
                 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)));