Vectorize tfer* functions

+ Harmonize common part of tfer*diff functions into a single tfer_diff function taking G0 as an input.
+ Add test.main_vec.
+ tfer_W1_pb, updated to use parse_inputs rs value.
+ Spacing updates within functions.

parse_inputs.m

@@ -21,7 +21,7 @@
 %           the integer charge state and particle mobility. 
 %=========================================================================%
 
-function [tau,C0,D,rs] = parse_inputs(sp,m,d,z,prop)
+function [tau, C0, D, rs] = parse_inputs(sp, m, d, z, prop)
 
 
 %-- Parse inputs ---------------------------------------------------------%
@@ -45,7 +45,7 @@
 %-- Evaluate output parameters -------------------------------------------%
 tau = B .* m;
 D = prop.D(B) .* z; % diffusion as a function of mechanical mobiltiy and charge state
-C0 = sp.V .* q ./ log(1/prop.r_hat); % calcualte recurring C0 parameter
+C0 = [sp.V]' .* q ./ log(1/prop.r_hat); % calcualte recurring C0 parameter
 
 if nargout>=4 % if required, calculate equilbirium radius
 % Note: Whether to pick the +ive of -ive root for rs is chosen based on a
@@ -54,22 +54,22 @@
 % conditions, where the +ive root should always be used).
 
     % evaluate +ive and -ive roots
-    r_m = (sqrt(C0./m) - ...
-        sqrt(C0./m - 4*sp.alpha*sp.beta)) ./ (2*sp.alpha);
-    r_p = (sqrt(C0./m) + ...
-        sqrt(C0./m - 4*sp.alpha*sp.beta)) ./ (2*sp.alpha);
+    r_m = (sqrt(C0 ./ m) - ...
+        sqrt(C0 ./ m - 4 .* [sp.alpha]' .* [sp.beta]')) ./ (2 .* [sp.alpha]');
+    r_p = (sqrt(C0 ./ m) + ...
+        sqrt(C0 ./ m - 4 .* [sp.alpha]' .* [sp.beta]')) ./ (2 .* [sp.alpha]');
 
     % determine which root is closer to centerline radius
-    [~,idx] = min([abs(r_m-prop.rc); abs(r_p-prop.rc)]);
-    idx(r_m==0) = 2; % avoid zero values for APM case
+    bo = abs(r_m-prop.rc) > abs(r_p-prop.rc);
+    bo(r_m==0) = 1; % avoid zero values for APM case
 
     % assign one of the two roots to rs
     rs = r_m; % by default use -ive root
-    rs(idx==2) = r_p(idx==2); % if closer to +ive root, use +ive root
+    rs(bo) = r_p(bo); % if closer to +ive root, use +ive root
 
 
     % zero out cases where no equilibrium radius exists (also removes complex numbers)
-    rs(C0./m < (4*sp.alpha*sp.beta)) = 0;
+    rs(C0 ./ m < (4 .* [sp.alpha]' .* [sp.beta]')) = 0;
 end
 
 

test/main_vec.m

@@ -0,0 +1,51 @@
+
+% MAIN_VEC  A function to test vector evaluation of transfer functions.
+%  AUTHOR: Timothy Sipkens, 2021-10-25
+
+%-- Initialize script ----------------------------------------------------%
+clear;
+close all;
+
+Rm = 3; % equivalent resolution of transfer functions (Reavell et al.)
+
+m_star = logspace(-2, 0, 7) .* 1e-18; % mass in kg (1 fg = 1e-18 kg)
+m = logspace(-2.3, 0.3, 801) .* 1e-18; % vector of mass
+
+z = 1; % integer charge state
+
+rho_eff100 = 900; % effective density
+Dm = 3;
+m100 = rho_eff100 * (pi * (100e-9)^3 / 6);
+m0 = m100 * (1/100) ^ Dm;
+d = 1e-9 .* (m ./ m0) .^ (1/Dm);
+    % specify mobility diameter vector with constant effective density
+
+prop = prop_pma('olfert'); % get properties of the CPMA
+prop.omega_hat = 1;
+prop.m0 = m0; % copy mass-mobility relation info (only used to find Rm)
+prop.Dm = Dm;
+
+% prop.omega_hat = 1; % NOTE: Uncomment for APM condition
+
+sp = get_setpoint(prop, 'm_star', m_star, 'Rm', Rm);
+    % get setpoint parameters
+
+
+K0 = tfer_1C(sp, m, d, z, prop);
+K1 = tfer_1C_diff(sp, m, d, z, prop);
+
+%{
+K2 = [];
+for ii=1:length(sp)
+    K2(ii,:) = tfer_FD(sp(ii), m, d, z, prop);
+end
+%}
+
+figure(1);
+semilogx(m, K0);
+hold on;
+semilogx(m, K1, 'k--');
+% semilogx(m, K2, 'r--');
+hold off;
+
+

tfer_1C.m

@@ -16,25 +16,26 @@
 %   G0          Function mapping final to initial radial position
 %=========================================================================%
 
-function [Lambda,G0] = tfer_1C(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_1C(sp, m, d, z, prop)
 
-[tau,C0] = parse_inputs(sp,m,d,z,prop);
+[tau, C0] = parse_inputs(sp, m, d, z, prop);
         % parse inputs for common parameters
 
 %-- Taylor series expansion constants ------------------------------------%
-C3 = tau .* (sp.alpha^2 * prop.rc + ...
-    2 * sp.alpha * sp.beta / prop.rc + sp.beta^2 / (prop.rc^3) - C0./(m.*prop.rc));
-C4 = tau .* (sp.alpha^2 - 2 * sp.alpha * sp.beta / (prop.rc^2) - 3*sp.beta^2/(prop.rc^4)+C0./(m.*(prop.rc^2)));
+C3 = tau .* ([sp.alpha]' .^ 2 * prop.rc + ...
+    2 .* [sp.alpha]' .* [sp.beta]' ./ prop.rc + [sp.beta]' .^ 2 ./ (prop.rc^3) - C0 ./ (m .* prop.rc));
+C4 = tau .* ([sp.alpha]' .^ 2 - 2 .* [sp.alpha]' .* [sp.beta]' ./ (prop.rc^2) - ...
+    3 .* [sp.beta]' .^ 2 / (prop.rc^4) + C0 ./ (m .* (prop.rc^2)));
 
 
 %-- Evaluate G0 and transfer function ------------------------------------%
-G0 = @(r) prop.rc + (r - prop.rc + C3./C4) .* ...
+G0 = @(r) prop.rc + (r - prop.rc + C3 ./ C4) .* ...
     exp(-C4 .* prop.L ./ prop.v_bar) - C3./C4;
 
-ra = min(prop.r2,max(prop.r1,G0(prop.r1)));
-rb = min(prop.r2,max(prop.r1,G0(prop.r2)));
+ra = min(prop.r2, max(prop.r1, G0(prop.r1)));
+rb = min(prop.r2, max(prop.r1, G0(prop.r2)));
 
-Lambda = (1/(2*prop.del)).*(rb-ra);
+Lambda = (1 ./ (2 .* prop.del)) .* (rb - ra);
 
 end
 

tfer_1C_diff.m

@@ -16,29 +16,14 @@
 %   G0          Function mapping final to initial radial position
 %=========================================================================%
 
-function [Lambda,G0] = tfer_1C_diff(sp,m,d,z,prop)
-
-[~,~,D] = parse_inputs(sp,m,d,z,prop); % get diffusion coeff.
-sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter
+function [Lambda, G0] = tfer_1C_diff(sp, m, d, z, prop)
 
+[~, ~, D] = parse_inputs(sp, m, d, z, prop); % get diffusion coeff.
 
 %-- Evaluate relevant functions ------------------------------------------%
-[~,G0] = tfer_1C(sp,m,d,z,prop);
+[~, G0] = tfer_1C(sp, m, d, z, prop);
     % get G0 function for this case
 
-rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % recurring quantity
-kap_fun = @(G,r) ...
-    (G-r).*erf(rho_fun(G,r))+...
-    sig.*sqrt(2/pi).*exp(-rho_fun(G,r).^2); % define function for kappa
-
-
-%-- Evaluate the transfer function and its terms -------------------------%
-K22 = kap_fun(G0(prop.r2),prop.r2);
-K21 = kap_fun(G0(prop.r2),prop.r1);
-K12 = kap_fun(G0(prop.r1),prop.r2);
-K11 = kap_fun(G0(prop.r1),prop.r1);
-Lambda = -1/(4*prop.del).*(K22-K12-K21+K11);
-Lambda(K22>1e2) = 0; % remove cases with large values out of error fun. eval.
-Lambda(abs(Lambda)<1e-10) = 0; % remove cases with roundoff error
+Lambda = tfer_diff(G0, D, prop);  % apply general diffusive form
 
 end

tfer_1C_pb.m

@@ -16,33 +16,44 @@
 %   G0          Function mapping final to initial radial position
 %=========================================================================%
 
-function [Lambda,G0] = tfer_1C_pb(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_1C_pb(sp, m, d, z, prop)
 
 [tau,C0,~,rs] = parse_inputs(sp,m,d,z,prop);
         % parse inputs for common parameters
 
 %-- Taylor series expansion constants ------------------------------------%
-C3 = tau.*(sp.alpha^2*prop.rc+2*sp.alpha*sp.beta/prop.rc+sp.beta^2/(prop.rc^3)-C0./(m.*prop.rc));
-C4 = tau.*(sp.alpha^2-2*sp.alpha*sp.beta/(prop.rc^2)-3*sp.beta^2/(prop.rc^4)+C0./(m.*(prop.rc^2)));
-
-A1 = -3*prop.v_bar./(4.*C4.^3.*prop.del^2);
-A2 = 2.*(C3.^2-C4.^2.*prop.del^2);
-A3 = @(r,ii) C4(ii).^2.*(r-prop.rc).^2-2.*C3(ii).*C4(ii).*(r-prop.rc);
-
-
-%-- Set up F function for minimization -----------------------------------%
-F = @(r,ii) A1(ii).*(A2(ii).*log(abs(C4(ii).*(r-prop.rc)+C3(ii)))+A3(r,ii));
-min_fun = @(rL,r0,ii) F(rL,ii)-F(r0,ii)-prop.L;
-
-
-%-- Evaluate G0 and transfer function ------------------------------------%
-G0 = @(r) G_fun(min_fun,r,rs,prop.r1,prop.r2,sp.alpha,sp.beta);
-
-ra = min(prop.r2,max(prop.r1,G0(prop.r1)));
-rb = min(prop.r2,max(prop.r1,G0(prop.r2)));
-
-Lambda = 3/4.*(rb-ra)./prop.del-1/4.*((rb-prop.rc)./prop.del).^3+...
-    1/4.*((ra-prop.rc)./prop.del).^3;
+C3 = tau .* ([sp.alpha]' .^ 2 .* prop.rc + ...
+    2 .* [sp.alpha]' .* [sp.beta]' ./ prop.rc + ...
+    [sp.beta]' .^ 2 ./ (prop.rc .^ 3) - C0 ./ (m .* prop.rc));
+C4 = tau .* ([sp.alpha]' .^ 2 - ...
+    2 .* [sp.alpha]' .* [sp.beta]' ./ (prop.rc .^ 2) - ...
+    3 .* [sp.beta]' .^ 2 ./ (prop.rc .^ 4) + ...
+    C0 ./ (m .* (prop.rc .^ 2)));
+
+A1 = -3 * prop.v_bar ./ (4 .* C4 .^ 3 .* prop.del ^ 2);
+A2 = 2 .* (C3 .^ 2 - C4 .^ 2 .* prop.del ^ 2);
+
+% Loop over setpoints.
+Lambda = zeros(size(A1));
+for jj=1:size(Lambda,1)
+    A3 = @(r,ii) C4(jj,ii) .^ 2 .* (r-prop.rc) .^ 2 - 2 .* C3(jj,ii) .* C4(jj,ii) .* (r-prop.rc);
+
+    %-- Set up F function for minimization -------------------------------%
+    F = @(r,ii) A1(jj,ii) .* (A2(jj,ii) .* ...
+        log(abs(C4(jj,ii) .* (r - prop.rc) + C3(jj,ii))) + A3(r,ii));
+    min_fun = @(rL,r0,ii) F(rL,ii) - F(r0,ii) - prop.L;
+
+    %-- Evaluate G0 and transfer function --------------------------------%
+    G0 = @(r) G_fun(min_fun, r, rs(jj,:), ...
+        prop.r1, prop.r2, sp(jj).alpha, sp(jj).beta);
+
+    ra = min(prop.r2, max(prop.r1, G0(prop.r1)));
+    rb = min(prop.r2, max(prop.r1, G0(prop.r2)));
+
+    Lambda(jj,:) = 3/4 .* (rb - ra) ./ prop.del - ...
+        1/4 .* ((rb - prop.rc) ./ prop.del) .^ 3 + ...
+        1/4 .* ((ra - prop.rc) ./ prop.del) .^ 3;
+end
 
 end
 

tfer_1S.m

@@ -16,13 +16,13 @@
 %   G0          Function mapping final to initial radial position
 %=========================================================================%
 
-function [Lambda,G0] = tfer_1S(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_1S(sp, m, d, z, prop)
 
-[tau,~,~,rs] = parse_inputs(sp, m, d, z, prop);
+[tau, ~, ~, rs] = parse_inputs(sp, m, d, z, prop);
         % parse inputs for common parameters
 
 %-- Estimate device parameter --------------------------------------------%
-lam = 2 .* tau .* (sp.alpha^2 - sp.beta^2 ./ (rs.^4)) .* prop.L ./ prop.v_bar;
+lam = 2 .* tau .* ([sp.alpha]' .^ 2 - [sp.beta]' .^ 2 ./ (rs .^ 4)) .* prop.L ./ prop.v_bar;
 
 
 %-- Evaluate G0 and transfer function ------------------------------------%

tfer_1S_diff.m

@@ -19,26 +19,11 @@
 function [Lambda,G0] = tfer_1S_diff(sp,m,d,z,prop)
 
 [~,~,D] = parse_inputs(sp,m,d,z,prop); % get diffusion coeff.
-sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter
-
 
 %-- Evaluate relevant functions ------------------------------------------%
 [~,G0] = tfer_1S(sp,m,d,z,prop);
     % get G0 function for this case
 
-rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % recurring quantity
-kap_fun = @(G,r) ...
-    (G-r).*erf(rho_fun(G,r))+...
-    sig.*sqrt(2/pi).*exp(-rho_fun(G,r).^2); % define function for kappa
-
-
-%-- Evaluate the transfer function and its terms -------------------------%
-K22 = kap_fun(G0(prop.r2),prop.r2);
-K21 = kap_fun(G0(prop.r2),prop.r1);
-K12 = kap_fun(G0(prop.r1),prop.r2);
-K11 = kap_fun(G0(prop.r1),prop.r1);
-Lambda = -1/(4*prop.del).*(K22-K12-K21+K11);
-Lambda(K22>1e2) = 0; % remove cases with large values out of error fun. eval.
-Lambda(abs(Lambda)<1e-10) = 0; % remove cases with roundoff error
+Lambda = tfer_diff(G0, D, prop);  % apply general diffusive form
 
 end

tfer_1S_pb.m

@@ -22,26 +22,30 @@
         % parse inputs for common parameters
 
 %-- Estimate device parameter --------------------------------------------%
-lam = 2 .* tau .* (sp.alpha^2 - sp.beta^2 ./ (rs.^4)) .* prop.L ./ prop.v_bar;
+lam = 2 .* tau .* ([sp.alpha]' .^ 2 - [sp.beta]' .^ 2 ./ (rs.^4)) .* prop.L ./ prop.v_bar;
 
 A1 = -3 * prop.L ./ (2 .* lam .* prop.del^2);
 A2 = rs.^2 + prop.rc^2 - 2*prop.rc.*rs - prop.del^2;
-A3 = @(r,ii) (r.^2) ./2 + (rs(ii) - 2 * prop.rc) .* r;
 
+% Loop over setpoints.
+Lambda = zeros(size(A1));
+for jj=1:size(Lambda,1)
+    A3 = @(r,ii) (r.^2) ./2 + (rs(jj,ii) - 2 * prop.rc) .* r;
+
+    %-- Set up F function for minimization -------------------------------%
+    F = @(r,ii) A1(jj,ii) .* (A2(jj,ii) .* log(r - rs(jj,ii)) + A3(r,ii));
+    min_fun = @(rL,r0,ii) F(rL,ii) - F(r0,ii) - prop.L;
 
-%-- Set up F function for minimization -----------------------------------%
-F = @(r,ii) A1(ii) .* (A2(ii) .* log(r - rs(ii)) + A3(r,ii));
-min_fun = @(rL,r0,ii) F(rL,ii) - F(r0,ii) - prop.L;
+    %-- Evaluate G0 and transfer function --------------------------------%
+    G0 = @(r) G_fun(min_fun, r, rs(jj,:), ...
+        prop.r1, prop.r2, sp(jj).alpha, sp(jj).beta);
 
+    ra = min(prop.r2, max(prop.r1, G0(prop.r1)));
+    rb = min(prop.r2, max(prop.r1, G0(prop.r2)));
 
-%-- Evaluate G0 and transfer function ------------------------------------%
-G0 = @(r) G_fun(min_fun, r, rs, prop.r1, prop.r2, sp.alpha, sp.beta);
-
-ra = min(prop.r2, max(prop.r1, G0(prop.r1)));
-rb = min(prop.r2, max(prop.r1, G0(prop.r2)));
-
-Lambda = 3/4 .* (rb - ra) ./ prop.del - 1/4 .* ((rb - prop.rc) ./ prop.del).^3+...
-    1/4 .* ((ra - prop.rc) ./ prop.del).^3;
+    Lambda(jj,:) = 3/4 .* (rb - ra) ./ prop.del - 1/4 .* ((rb - prop.rc) ./ prop.del).^3+...
+        1/4 .* ((ra - prop.rc) ./ prop.del).^3;
+end
 
 end
 

tfer_2C.m

@@ -16,29 +16,36 @@
 %   G0          Function mapping final to initial radial position
 %=========================================================================%
 
-function [Lambda,G0] = tfer_2C(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_2C(sp, m, d, z, prop)
 
 [tau,C0] = parse_inputs(sp,m,d,z,prop);
         % parse inputs for common parameters
 
 %-- Taylor series expansion constants ------------------------------------%
-C3 = tau.*(sp.alpha^2*prop.rc+2*sp.alpha*sp.beta/prop.rc+sp.beta^2/(prop.rc^3)-C0./(m.*prop.rc));
-C4 = tau.*(sp.alpha^2-2*sp.alpha*sp.beta/(prop.rc^2)-3*sp.beta^2/(prop.rc^4)+C0./(m.*(prop.rc^2)));
-C5 = 2.*tau.*(2*sp.alpha*sp.beta./(prop.rc^3)+6*sp.beta^2/(prop.rc^5)-C0./(m.*(prop.rc^3)));
+C3 = tau.*([sp.alpha]' .^ 2 .* prop.rc + ...
+    2 .* [sp.alpha]' .* [sp.beta]' ./ prop.rc + ...
+    [sp.beta]' .^ 2 ./ (prop.rc .^ 3) - ...
+    C0 ./ (m .* prop.rc));
+C4 = tau .* ([sp.alpha]' .^ 2 - ...
+    2 .* [sp.alpha]' .* [sp.beta]' ./ (prop.rc .^ 2) - ...
+    3 .* [sp.beta]' .^ 2 ./ (prop.rc .^ 4) + C0 ./ (m .* (prop.rc .^ 2)));
+C5 = 2 .* tau .* (2 .* [sp.alpha]' .* [sp.beta]' ./ (prop.rc .^ 3) + ...
+    6 .* [sp.beta]' .^ 2 ./ (prop.rc .^ 5) - C0 ./ (m .* (prop.rc .^ 3)));
 
-C6 = -1./(sqrt(4.*C3.*C5-C4.^2));
+C6 = -1 ./ (sqrt(4 .* C3 .* C5 - C4 .^ 2));
 
-f = @(r) 2.*C6.*prop.v_bar.*...
-    atan(C6.*(2.*C5.*(r-prop.rc)+C4));
+f = @(r) 2 .* C6 .* prop.v_bar .* ...
+    atan(C6 .* (2 .* C5 .* (r - prop.rc) + C4));
 
 
 %-- Evaluate G0 and transfer function ------------------------------------%
-G0 = @(r) (tan((f(r)-prop.L)./(2.*C6.*prop.v_bar))./C6-C4)./(2.*C5)+prop.rc;
+G0 = @(r) (tan((f(r) - prop.L) ./ ...
+    (2 .* C6 .* prop.v_bar)) ./ C6 - C4) ./ (2.*C5)+prop.rc;
 
-ra = min(prop.r2,max(prop.r1,G0(prop.r1)));
-rb = min(prop.r2,max(prop.r1,G0(prop.r2)));
+ra = min(prop.r2, max(prop.r1, G0(prop.r1)));
+rb = min(prop.r2, max(prop.r1, G0(prop.r2)));
 
-Lambda = real((1/(2*prop.del)).*(rb-ra));
+Lambda = real((1 ./ (2 .* prop.del)) .* (rb - ra));
 
 end
 

tfer_2C_diff.m

@@ -26,19 +26,6 @@
 [~,G0] = tfer_2C(sp,m,d,z,prop);
     % get G0 function for this case
 
-rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % recurring quantity
-kap_fun = @(G,r) ...
-    (G-r).*erf(rho_fun(G,r))+...
-    sig.*sqrt(2/pi).*exp(-rho_fun(G,r).^2); % define function for kappa
-
-
-%-- Evaluate the transfer function and its terms -------------------------%
-K22 = kap_fun(real(G0(prop.r2)),prop.r2);
-K21 = kap_fun(real(G0(prop.r2)),prop.r1);
-K12 = kap_fun(real(G0(prop.r1)),prop.r2);
-K11 = kap_fun(real(G0(prop.r1)),prop.r1);
-Lambda = -1/(4*prop.del).*(K22-K12-K21+K11);
-Lambda(K22>1e2) = 0; % remove cases with large values out of error fun. eval.
-Lambda(abs(Lambda)<1e-10) = 0; % remove cases with roundoff error
+Lambda = tfer_diff(G0, D, prop);  % apply general diffusive form
 
 end