Refactors Abeles Calculation

targetFunctions/common/reflectivityCalculations/abeles/abelesParallelPoints.m

@@ -1,84 +1,55 @@
 function ref = abelesParallelPoints(q,N,layersThick,layersRho,layersSigma)
 
-tiny = 1e-30;
-ci = complex(0,1);
-c0 = complex(0,0);
+% Vectorised version of reflectivity with complex rho
 
+c1 = complex(1,0);
 ref = zeros(length(q),1);
 
-layersRho1 = layersRho(1);
-layersRho2 = layersRho(2);
-layersSigma2 = layersSigma(2);
+% SLD is always relative to bulk in value
+sld = layersRho(2:end) - layersRho(1);
+
+% Always need the next sigma squared value for layers 1:end-1
+nextSigmaSquared = layersSigma(2:end).^2;
 
 parfor points = 1:length(q)
 
-    M_tot = [c0 c0 ; c0 c0];
-    M_n = [c0 c0 ; c0 c0];
-    M_res = [c0 c0 ; c0 c0];
-
-    Q = q(points);
-
-    bulkInSLD = layersRho1;
-    bulkInSLD = bulkInSLD + complex(0,tiny);
-
-    k0 = 0.5 * Q;
-
-    % Find k1..
-    sld_1 = layersRho2 - bulkInSLD;
-    k1 = findkn(k0, sld_1);
-
-    % Find r01
-    nom1 = k0 - k1;
-    denom1 = k0 + k1;
-    sigmasqrd = layersSigma2 ^ 2;
-    err1 = exp(-2 * k1 * k0 * sigmasqrd);
-    r01 = (nom1 / denom1) * err1;
-
-    % Generate the M1 matrix:
-    M_tot(1,1) = complex(1,0);
-    M_tot(1,2) = r01;
-    M_tot(2,1) = r01;
-    M_tot(2,2) = complex(1,0);
-
-    kn_ptr = k1;
-
-    for n = 2:N-1
-
-        % Find kn and k_n+1 (ex. k1 and k2 for n=1): */
-        sld_np1 = layersRho(n + 1);
-        sld_np1 = sld_np1 - bulkInSLD;
-
-        kn = kn_ptr;
-        knp1 = findkn(k0, sld_np1);
-
-        % Find r_n,n+1:
-        nom_n = kn - knp1;
-        denom_n = kn + knp1;
-        sigmasqrd = layersSigma(n + 1)^2;
-        err_n = exp(-2 * kn * knp1 * sigmasqrd);
-        r_n_np1 = (nom_n / denom_n) * err_n;
-
+    M = eye(2, 'like', c1);
+    k0 = 0.5 * q(points);
+    kn = complex(k0);
+
+    for n = 1:N-1
+
+        knp1 = findkn(k0, sld(n));
+
         % Find the Phase Factor = (k_n * d_n)
-        beta = kn * layersThick(n) * ci;
-
-        % Create the M_n matrix: */
-        M_n(1,1) = exp(beta);
-        M_n(1,2) = r_n_np1 * exp(beta);
-        M_n(2,1) = r_n_np1 * exp(-beta);
-        M_n(2,2) = exp(-beta);
-
-        % Multiply the matrices
-        M_res = M_tot * M_n;
-
-        % Reassign the values back to M_tot:
-        M_tot = M_res;
-
-        % Point to k_n+1 and sld_n+1 via kn_ptr sld_n_ptr:
-        kn_ptr = knp1;
+        beta = kn * layersThick(n);
+
+        numerator = knp1 - kn;
+        denominator = knp1 + kn;
+        erf = exp(-2 * knp1 * kn * nextSigmaSquared(n));
+        r_n = (numerator / denominator) * erf;
+
+        % Multiply system transfer matrix by the layer transfer matrix
+        % Coder behaves better if you just do it manually
+        % rather than M = M * M_n;
+        % we'll denote the entries M = [a b; c d]
+        exp_beta = exp(1i*beta);
+        inv_exp_beta = exp(-1i*beta);
+
+        a_exp_beta = M(1,1) * exp_beta;
+        b_inv_exp_beta = M(1,2) * inv_exp_beta;
+        c_exp_beta = M(2,1) * exp_beta;
+        d_inv_exp_beta = M(2,2) * inv_exp_beta;
+
+        M = [a_exp_beta + r_n*b_inv_exp_beta r_n*a_exp_beta + b_inv_exp_beta;
+             c_exp_beta + r_n*d_inv_exp_beta r_n*c_exp_beta + d_inv_exp_beta];
+
+        % Set value of kn for the next layer
+        kn = knp1;
 
     end
 
-    R = abs(M_res(2,1) / M_res(1,1));
+    R = abs(M(2,1) / M(1,1));
     ref(points) = R^2;
 
 end

targetFunctions/common/reflectivityCalculations/abeles/abelesSingle.m

@@ -1,82 +1,55 @@
 function ref = abelesSingle(q,N,layersThick,layersRho,layersSigma)
 
-% New Matlab version of reflectivity
-% with complex rho...
-
-% Pre-allocation
-tiny = 1e-30;
-ci = complex(0,1);
-c0 = complex(0,0);
-M_tot = [c0 c0 ; c0 c0];
-M_n = [c0 c0 ; c0 c0];
-M_res = [c0 c0 ; c0 c0];
-ref = zeros(length(q),1);
-
-for points = 1:length(q)
-
-    Q = q(points);
-
-    bulkInSLD = layersRho(1);
-    bulkInSLD = bulkInSLD + complex(0,tiny);
-
-    k0 = 0.5 * Q;
-
-    % Find k1..
-    sld_1 = layersRho(2) - bulkInSLD;
-    k1 = findkn(k0, sld_1);
+% Vectorised version of reflectivity with complex rho
 
-    % Find r01
-    nom1 = k0 - k1;
-    denom1 = k0 + k1;
-    sigmasqrd = layersSigma(2) ^ 2;
-    err1 = exp(-2 * k1 * k0 * sigmasqrd);
-    r01 = (nom1 / denom1) * err1;
+c1 = complex(1,0);
+ref = zeros(length(q),1);
 
-    % Generate the M1 matrix:
-    M_tot(1,1) = complex(1,0);
-    M_tot(1,2) = r01;
-    M_tot(2,1) = r01;
-    M_tot(2,2) = complex(1,0);
+% SLD is always relative to bulk in value
+sld = layersRho(2:end) - layersRho(1);
 
-    kn_ptr = k1;
+% Always need the next sigma squared value for layers 1:end-1
+nextSigmaSquared = layersSigma(2:end).^2;
 
-    for n = 2:N-1
+for points = 1:length(q)
 
-        % Find kn and k_n+1 (ex. k1 and k2 for n=1): */
-        sld_np1 = layersRho(n + 1);
-        sld_np1 = sld_np1 - bulkInSLD;
+    M = eye(2, 'like', c1);
+    k0 = 0.5 * q(points);
+    kn = complex(k0);
 
-        kn = kn_ptr;
-        knp1 = findkn(k0, sld_np1);
+    for n = 1:N-1
 
-        % Find r_n,n+1:
-        nom_n = kn - knp1;
-        denom_n = kn + knp1;
-        sigmasqrd = layersSigma(n + 1)^2;
-        err_n = exp(-2 * kn * knp1 * sigmasqrd);
-        r_n_np1 = (nom_n / denom_n) * err_n;
+        knp1 = findkn(k0, sld(n));
 
         % Find the Phase Factor = (k_n * d_n)
-        beta = kn * layersThick(n) * ci;
-
-        % Create the M_n matrix: */
-        M_n(1,1) = exp(beta);
-        M_n(1,2) = r_n_np1 * exp(beta);
-        M_n(2,1) = r_n_np1 * exp(-beta);
-        M_n(2,2) = exp(-beta);
-
-        % Multiply the matrices
-        M_res = M_tot * M_n;
-
-        % Reassign the values back to M_tot:
-        M_tot = M_res;
-
-        % Point to k_n+1 and sld_n+1 via kn_ptr sld_n_ptr:
-        kn_ptr = knp1;
+        beta = kn * layersThick(n);
+
+        numerator = knp1 - kn;
+        denominator = knp1 + kn;
+        erf = exp(-2 * knp1 * kn * nextSigmaSquared(n));
+        r_n = (numerator / denominator) * erf;
+
+        % Multiply system transfer matrix by the layer transfer matrix
+        % Coder behaves better if you just do it manually
+        % rather than M = M * M_n;
+        % we'll denote the entries M = [a b; c d]
+        exp_beta = exp(1i*beta);
+        inv_exp_beta = exp(-1i*beta);
+
+        a_exp_beta = M(1,1) * exp_beta;
+        b_inv_exp_beta = M(1,2) * inv_exp_beta;
+        c_exp_beta = M(2,1) * exp_beta;
+        d_inv_exp_beta = M(2,2) * inv_exp_beta;
+
+        M = [a_exp_beta + r_n*b_inv_exp_beta r_n*a_exp_beta + b_inv_exp_beta;
+             c_exp_beta + r_n*d_inv_exp_beta r_n*c_exp_beta + d_inv_exp_beta];
+
+        % Set value of kn for the next layer
+        kn = knp1;
 
     end
 
-    R = abs(M_res(2,1) / M_res(1,1));
+    R = abs(M(2,1) / M(1,1));
     ref(points) = R^2;
 
 end

targetFunctions/common/resolutionFunctions/simpleGaussian/resolutionPolly.m

@@ -14,10 +14,8 @@
         sumg = sumg + g;
         simulation(j) = simulation(j) + rawSimulation(i+j) * g;
     end
-
-    if (sumg ~= 0)
-        simulation(j) = simulation(j) / sumg;
-    end
+
+    simulation(j) = simulation(j) / sumg;
 
 end