diff --git a/parse_inputs.m b/parse_inputs.m
index 7889b91..4f3fb63 100644
--- a/parse_inputs.m
+++ b/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
diff --git a/test/main_vec.m b/test/main_vec.m
new file mode 100644
index 0000000..1cf8d32
--- /dev/null
+++ b/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;
+
+
diff --git a/tfer_1C.m b/tfer_1C.m
index 29b1dbe..3959f63 100644
--- a/tfer_1C.m
+++ b/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
diff --git a/tfer_1C_diff.m b/tfer_1C_diff.m
index c2043ef..2f869c0 100644
--- a/tfer_1C_diff.m
+++ b/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
diff --git a/tfer_1C_pb.m b/tfer_1C_pb.m
index 5915564..7c14123 100644
--- a/tfer_1C_pb.m
+++ b/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
diff --git a/tfer_1S.m b/tfer_1S.m
index fca4bc1..552be47 100644
--- a/tfer_1S.m
+++ b/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 ------------------------------------%
diff --git a/tfer_1S_diff.m b/tfer_1S_diff.m
index 2dd6d10..369871a 100644
--- a/tfer_1S_diff.m
+++ b/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
diff --git a/tfer_1S_pb.m b/tfer_1S_pb.m
index f8a9b0a..44b3c1c 100644
--- a/tfer_1S_pb.m
+++ b/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
diff --git a/tfer_2C.m b/tfer_2C.m
index f41ee33..49ee233 100644
--- a/tfer_2C.m
+++ b/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
diff --git a/tfer_2C_diff.m b/tfer_2C_diff.m
index 34e7cd7..4ba8d1f 100644
--- a/tfer_2C_diff.m
+++ b/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
diff --git a/tfer_2S.m b/tfer_2S.m
index 747475b..071b1f7 100644
--- a/tfer_2S.m
+++ b/tfer_2S.m
@@ -16,28 +16,31 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_2S(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_2S(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;
%-- Taylor series expansion constants ------------------------------------%
-C1 = 2.*tau.*(sp.alpha^2-sp.beta^2./(rs.^4));
-C2 = -2.*tau.*(sp.alpha^2./rs-5*sp.beta^2./(rs.^5));
+C1 = 2 .* tau .* ([sp.alpha]' .^ 2 - ...
+ [sp.beta]' .^ 2 ./ (rs .^ 4));
+C2 = -2 .* tau .* ([sp.alpha]' .^ 2 ./ rs - ...
+ 5 .* [sp.beta]' .^ 2 ./ (rs .^ 5));
%-- Evaluate G0 and transfer function ------------------------------------%
-G0 = @(r) rs+C1.*(r-rs).*exp(-lam)./...
- (C2.*(r-rs)+C1-C2.*(r-rs).*exp(-lam));
+G0 = @(r) rs + C1 .* (r - rs) .* exp(-lam) ./ ...
+ (C2 .* (r - rs) + C1 - C2 .* (r - rs) .* exp(-lam));
-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
diff --git a/tfer_2S_diff.m b/tfer_2S_diff.m
index 5c4d6ef..f833d24 100644
--- a/tfer_2S_diff.m
+++ b/tfer_2S_diff.m
@@ -19,26 +19,11 @@
function [Lambda,G0] = tfer_2S_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_2S(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
diff --git a/tfer_FD.m b/tfer_FD.m
index 0927980..4da2ddc 100644
--- a/tfer_FD.m
+++ b/tfer_FD.m
@@ -22,101 +22,102 @@
function [Lambda,n] = tfer_FD(sp,m,d,z,prop)
-[tau,C0,D] = parse_inputs(sp,m,d,z,prop);
+[tau, C0, D] = parse_inputs(sp,m,d,z,prop);
% parse inputs for common parameters
%-- Discretize the space -------------------------------------------------%
nr = 200;
-dr = (prop.r2-prop.r1)/nr;
-r_vec = (prop.r1+dr/2):dr:(prop.r2-dr/2);
+dr = (prop.r2 - prop.r1) / nr;
+r_vec = (prop.r1 + dr/2):dr:(prop.r2 - dr/2);
nz = 200;
-dz = prop.L/(nz-1);
-if nargout>=2; n.z_vec=0:dz:prop.L; n.r_vec=r_vec; end % vector of z positions, used for plotting only
+dz = prop.L /( nz-1);
+if nargout >= 2; n.z_vec = 0:dz:prop.L; n.r_vec = r_vec; end % vector of z positions, used for plotting only
%-- Parameters related to CPMA geometry ----------------------------------%
-del = (prop.r2-prop.r1)/2;
-rc = (prop.r2+prop.r1)/2;
-
-D0 = D.*prop.L/(prop.del^2*prop.v_bar);
+D0 = D .* prop.L / (prop.del ^ 2 * prop.v_bar);
% dimensionless diffusion coefficient
%-- Evaluate relevant radial quantities ----------------------------------%
-v_theta = sp.alpha.*r_vec+sp.beta./r_vec; % azimuthal velocity distribution
-F_e = C0./r_vec; % electrostatic force
-v_z = 3/2*prop.v_bar.*(1-((r_vec-rc)./(del)).^2);
+v_theta = [sp.alpha]' .* r_vec + [sp.beta]' ./ r_vec; % azimuthal velocity distribution
+F_e = C0 ./ r_vec; % electrostatic force
+v_z = 3/2 .* prop.v_bar .* (1 - ((r_vec - prop.rc) ./ (prop.del)) .^ 2);
% axial velocity distribution (parabolic)
% v_z = prop.v_bar.*ones(size(r_vec)); % axial velocity distribution (plug)
-%-- Speed computation using resolution to limit computation --------------%
-ind = 1:length(m);
-if isfield(sp,'Rm') % if resolution is specified, use to reduce necessary computation
- cond0 = or(m>(z.*sp.m_star+2.*sp.m_max),...
- m<(z.*sp.m_star-2.*sp.m_max));
- % NOTE: conditions limits consideration of those regions where particles
- % do not escape, speeding computation.
- ind = ind(~cond0);
-end
+Lambda = zeros(length(sp), length(m));% initialize the transfer function variable
+for jj=1:length(sp)
+ %-- Speed computation using resolution to limit computation --------------%
+ ind = 1:length(m);
+ if isfield(sp(jj), 'Rm') % if resolution is specified, use to reduce necessary computation
+ cond0 = or(m > (z .* sp(jj).m_star + 2 .* sp(jj).m_max),...
+ m < (z .* sp(jj).m_star - 2 .* sp(jj).m_max));
+ % NOTE: conditions limits consideration of those regions where particles
+ % do not escape, speeding computation.
+ ind = ind(~cond0);
+ end
-%-- Loop over mass, but not m_star ---------------------------------------%
-Lambda = zeros(1,length(m));% initialize the transfer function variable
-for ii=ind % loop over mass (not m_star)
-
- F_c = m(ii).*v_theta.^2./r_vec; % centriputal force
- c_r = tau(ii)/m(ii).*(F_c-F_e); % particle velocity across streamlines
- drcr_dr = tau(ii).*(2*sp.alpha^2.*r_vec-2*sp.beta^2./(r_vec.^3));
+ %-- Loop over mass, but not m_star -----------------------------------%
+ for ii=ind % loop over mass (not m_star)
-
- %-- Initialize variables -----------------------------------%
- n_vec = ones(size(r_vec)); % number concentration at current axial position
- n_vec0 = n_vec; % initial number concentration (used for func. eval.)
- if sp.m_star>=2 % initilize variables used for visualizing number concentrations
- n_mat = zeros(nz,length(r_vec));
- n_mat(1,:) = n_vec;
- end
-
-
- %-- Get coefficients -------------------------------------%
- zet = v_z./dz;
- gam = D(ii)/(2*dr^2);
- kap = D(ii)./(4.*r_vec.*dr);
- eta = 1./(2.*r_vec).*drcr_dr;
- phi = c_r./(4*dr);
-
- %== Crank-Nicholson ========================================%
- [a,b,c,RHS] = gen_eqs(zet,gam,kap,eta,phi);
-
- %-- Primary loop for finite difference ---------------------%
- for jj = 2:nz
- n_vec = tridiag([0,a],b,c,RHS(n_vec));
- % solve system using Thomas algorithm
+ F_c = m(ii) .* v_theta(jj,:) .^ 2 ./ r_vec; % centriputal force
+ c_r = tau(ii) / m(ii) .* (F_c - F_e(jj,:)); % particle velocity across streamlines
+ drcr_dr = tau(ii) .* (2 .* sp(jj).alpha .^ 2 .* r_vec - ...
+ 2 .* sp(jj).beta .^ 2 ./ (r_vec .^ 3));
- if D0(ii)<1e-3 % for low diffusion, stabalize by smoothing oscillations
- n_vec = [0,n_vec,0];
- n_vec = (n_vec(1:end-2)+1e2.*n_vec(2:end-1)+n_vec(3:end))./...
- (1e2+2);
- end
- if nargout>=2; n_mat(jj,:) = n_vec; end
- % record particle distribution at this axial position
- end
-
- %-- Format output ------------------------------------------%
- if nargout>=2 % for visualizing number concentrations
- n.n_mat{ii} = max(n_mat,0);
+ %-- Initialize variables -----------------------------------%
+ n_vec = ones(size(r_vec)); % number concentration at current axial position
+ n_vec0 = n_vec; % initial number concentration (used for func. eval.)
+ if sp(jj).m_star >= 2 % initilize variables used for visualizing number concentrations
+ n_mat = zeros(nz, length(r_vec));
+ n_mat(1,:) = n_vec;
+ end
+
+
+ %-- Get coefficients -------------------------------------%
+ zet = v_z ./ dz;
+ gam = D(ii) / (2 * dr ^ 2);
+ kap = D(ii) ./ (4 .* r_vec .* dr);
+ eta = 1 ./ (2 .* r_vec) .* drcr_dr;
+ phi = c_r ./ (4 * dr);
+
+ %== Crank-Nicholson ========================================%
+ [a, b, c, RHS] = gen_eqs(zet, gam, kap, eta, phi);
+
+ %-- Primary loop for finite difference ---------------------%
+ for kk = 2:nz
+ n_vec = tridiag([0,a], b, c, RHS(n_vec));
+ % solve system using Thomas algorithm
+
+ if D0(ii)<1e-3 % for low diffusion, stabalize by smoothing oscillations
+ n_vec = [0,n_vec,0];
+ n_vec = (n_vec(1:end-2) + 1e2 .* n_vec(2:end-1) + ...
+ n_vec(3:end)) ./ ...
+ (1e2 + 2);
+ end
+
+ if nargout>=2; n_mat(kk,:) = n_vec; end
+ % record particle distribution at this axial position
+ end
+
+ %-- Format output ------------------------------------------%
+ if nargout>=2 % for visualizing number concentrations
+ n.n_mat{jj,ii} = max(n_mat,0);
+ end
+
+ Lambda(jj,ii) = sum(n_vec .* v_z) ./ sum(n_vec0 .* v_z);
+ % evaluate transfer fucntion
+
+ if Lambda(jj,ii) < 0.0005; Lambda(jj,ii) = 0; end
+ % truncate small func. values
+
end
-
- Lambda(ii) = sum(n_vec.*v_z)/sum(n_vec0.*v_z);
- % evaluate transfer fucntion
-
- if Lambda(ii) < 0.0005; Lambda(ii) = 0; end
- % truncate small func. values
-
end
end
@@ -125,37 +126,36 @@
%== GEN_EQS ==============================================================%
% Generate matrix equations for Crank-Nicholson scheme
% Author: Timothy Sipkens, 2018-12-27
-function [a,b,c,RHS] = gen_eqs(zet,gam,kap,eta,phi)
+function [a, b, c, RHS] = gen_eqs(zet, gam, kap, eta, phi)
ind_mid = 2:(length(zet)-1);
%-- Righ-hand side of eq. to be solved --------------------%
-RHS1 = @(n) (zet(1)-2*gam-eta(1)).*n(1)+...
- (gam+kap(1)-phi(1)).*n(2);
-RHS2 = @(n) (zet(ind_mid)-2*gam-eta(ind_mid)).*n(2:(end-1))+...
- (gam-kap(ind_mid)+phi(ind_mid)).*n(1:(end-2))+...
- (gam+kap(ind_mid)-phi(ind_mid)).*n(3:end);
-RHS3 = @(n) (zet(end)-2*gam-eta(end)).*n(end)+...
- (gam-kap(end)+phi(end)).*n(end-1);
-RHS = @(n) [RHS1(n),RHS2(n),RHS3(n)];
+RHS1 = @(n) (zet(1) - 2 * gam - eta(1)).*n(1)+...
+ (gam+kap(1) - phi(1)) .* n(2);
+RHS2 = @(n) (zet(ind_mid) - 2 * gam-eta(ind_mid)) .* n(2:(end-1)) + ...
+ (gam - kap(ind_mid) + phi(ind_mid)) .* n(1:(end-2)) + ...
+ (gam + kap(ind_mid) - phi(ind_mid)) .* n(3:end);
+RHS3 = @(n) (zet(end) - 2 * gam - eta(end)) .* n(end)+...
+ (gam - kap(end) + phi(end)) .* n(end-1);
+RHS = @(n) [RHS1(n), RHS2(n), RHS3(n)];
%-- Form tridiagonal matrix ------------------------------%
-b1 = zet(1)+2*gam+eta(1); % center
-c1 = -gam-kap(1)+phi(1); % +1
-
-a2 = (-gam+kap(ind_mid)-phi(ind_mid)); % -1
-b2 = zet(ind_mid)+2*gam+eta(ind_mid); % center
-c2 = -gam-kap(ind_mid)+phi(ind_mid); % +1
+b1 = zet(1) + 2 .* gam + eta(1); % center
+c1 = -gam - kap(1) + phi(1); % +1
-a3 = -gam+kap(end)-phi(end); % -1
-b3 = zet(end)+2*gam+eta(end); % center
+a2 = (-gam + kap(ind_mid) - phi(ind_mid)); % -1
+b2 = zet(ind_mid) + 2 * gam + eta(ind_mid); % center
+c2 = -gam - kap(ind_mid) + phi(ind_mid); % +1
-a = [a2,a3];
-b = [b1,b2,b3];
-c = [c1,c2];
+a3 = -gam + kap(end) - phi(end); % -1
+b3 = zet(end) + 2 * gam+eta(end); % center
+a = [a2, a3];
+b = [b1, b2, b3];
+c = [c1, c2];
end
diff --git a/tfer_GE.m b/tfer_GE.m
index 1eaddb7..9729d32 100644
--- a/tfer_GE.m
+++ b/tfer_GE.m
@@ -16,32 +16,38 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_GE(sp,m,d,z,prop)
+function [Lambda,G0] = tfer_GE(sp, m, d, z, prop)
-[tau,C0,~,rs] = parse_inputs(sp,m,d,z,prop);
+[tau, C0, ~, rs] = parse_inputs(sp, m, d, z, prop);
% parse inputs for common parameters
%-- Estimate recurring quantities ----------------------------------------%
-C6 = 2*sp.alpha*sp.beta-C0./m;
-C7 = sqrt(4*sp.alpha^2*sp.beta^2-C6.^2);
+C6 = 2 .* [sp.alpha]' .* [sp.beta]' - C0 ./ m;
+C7 = sqrt(4 .* [sp.alpha]' .^ 2 .* [sp.beta]' .^ 2 - C6 .^ 2);
-A1 = prop.v_bar./(4.*tau.*sp.alpha^2);
-A2 = 2.*C6./C7;
+A1 = prop.v_bar ./ (4 .* tau .* [sp.alpha]' .^ 2);
+A2 = 2 .* C6 ./ C7;
-%-- Set up F function for minimization -----------------------------------%
-F = @(r,ii) A1(ii).*(log(sp.alpha^2.*r^4+sp.beta^2+C6(ii).*r.^2)-...
- A2(ii).*atan((2*sp.alpha^2.*r.^2+C6(ii))./C7(ii)));
-min_fun = @(rL,r0,ii) F(rL,ii)-F(r0,ii)-prop.L;
+Lambda = zeros(size(A1));
+for jj=1:size(Lambda,1)
+ %-- Set up F function for minimization -------------------------------%
+ F = @(r,ii) A1(jj,ii) .* (log(sp(jj).alpha .^ 2 .* r .^ 4 + ...
+ sp(jj).beta .^ 2 + C6(jj,ii) .* r .^ 2) - ...
+ A2(jj,ii) .* atan((2 .* sp(jj).alpha .^ 2 .* r .^ 2 + ...
+ C6(jj,ii)) ./ C7(jj,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);
+ %-- 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)));
-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(jj,:) = (1 / (2 .* prop.del)).*(rb - ra);
+end
end
diff --git a/tfer_GE_diff.m b/tfer_GE_diff.m
index aaf45e1..6b5435d 100644
--- a/tfer_GE_diff.m
+++ b/tfer_GE_diff.m
@@ -16,29 +16,17 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_GE_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_GE_diff(sp, m, d, z, prop)
+[~, ~, D] = parse_inputs(sp, m, d, z, prop); % get diffusion coeff.
%-- Evaluate relevant functions ------------------------------------------%
-[~,G0] = tfer_GE(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 = zeros(length(sp), length(m));
+for jj=1:size(Lambda,1)
+ [~, G0] = tfer_GE(sp(jj), m, d, z, prop);
+ % get G0 function for this case
+
+ Lambda(jj,:) = tfer_diff(G0, D, prop); % apply general diffusive form
+end
end
diff --git a/tfer_W1.m b/tfer_W1.m
index a2ce4b6..42a90a3 100644
--- a/tfer_W1.m
+++ b/tfer_W1.m
@@ -16,23 +16,29 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_W1(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_W1(sp, m, d, z, prop)
-[tau,C0,~,rs] = parse_inputs(sp,m,d,z,prop);
+if prop.omega_hat ~= 1
+ warning(['WARNING: omega_hat ~= 1! ', ...
+ 'tfer_W1 output is likely inaccurate.'])
+end
+
+[tau, C0, ~, 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 ------------------------------------%
-G0 = @(r) 1./(sp.omega1.*sqrt(m)).*...
- sqrt((m.*sp.omega1^2.*r.^2-C0).*exp(-lam)+C0);
+G0 = @(r) 1 ./ ([sp.omega1]' .* sqrt(m)).*...
+ sqrt((m .* [sp.omega1]' .^ 2 .* r .^ 2 - C0) .* exp(-lam) + C0);
-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
diff --git a/tfer_W1_diff.m b/tfer_W1_diff.m
index 68fc492..65d2bac 100644
--- a/tfer_W1_diff.m
+++ b/tfer_W1_diff.m
@@ -16,28 +16,14 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_W1_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_W1_diff(sp, m, d, z, prop)
+[~,~,D] = parse_inputs(sp, m, d, z, prop); % get diffusion coeff.
%-- Evaluate relevant functions ------------------------------------------%
-[~,G0] = tfer_W1(sp,m,d,z,prop);
+[~, G0] = tfer_W1(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
diff --git a/tfer_W1_pb.m b/tfer_W1_pb.m
index 2c434b4..89622ea 100644
--- a/tfer_W1_pb.m
+++ b/tfer_W1_pb.m
@@ -16,40 +16,41 @@
% G0 Function mapping final to initial radial position
%=========================================================================%
-function [Lambda,G0] = tfer_W1_pb(sp,m,d,z,prop)
+function [Lambda, G0] = tfer_W1_pb(sp, m, d, z, prop)
-[tau,C0] = parse_inputs(sp,m,d,z,prop);
+[tau, C0, ~, rs] = parse_inputs(sp, m, d, z, prop);
% parse inputs for common parameters
-%-- Estimate equilibrium radius ------------------------------------------%
-if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc
- rs = real((sqrt(C0./m)-sqrt(C0./m-4*sp.alpha*sp.beta))./(2*sp.alpha)); % equiblirium radius for a given mass
-else
- rs = real((sqrt(C0./m)+sqrt(C0./m-4*sp.alpha*sp.beta))./(2*sp.alpha)); % equiblirium radius for a given mass
-end
-
-
%-- Estimate recurring quantities ----------------------------------------%
-A1 = -3.*prop.v_bar./(4.*tau.*sp.omega1.^4.*prop.del^2);
-A2 = sp.omega1.^2.*(prop.rc^2-prop.del^2)+C0./m;
-A3 = 4*prop.rc.*sp.omega1.*sqrt(C0./m);
-A4 = @(r,ii) sp.omega1.^2.*(r.^2-4*prop.rc.*r);
-
-
-%-- Set up F function for minimization -----------------------------------%
-F = @(r,ii) A1(ii).*(A2(ii).*log(C0./m(ii)-(sp.omega1.*r).^2)+...
- A3(ii).*atanh(sqrt(m(ii)./C0).*sp.omega1.*r)+A4(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;
+A1 = -3 .* prop.v_bar ./ (4 .* tau .* [sp.omega1]' .^ 4 .* prop.del .^ 2);
+A2 = [sp.omega1]' .^ 2 .* (prop.rc .^ 2 - prop.del .^ 2) + C0 ./ m;
+A3 = 4 .* prop.rc .* [sp.omega1]' .* sqrt(C0 ./ m);
+
+
+% Loop over setpoints.
+Lambda = zeros(size(A1));
+for jj=1:size(Lambda,1)
+
+ A4 = @(r,ii) sp(jj).omega1 .^ 2 .* (r .^ 2 - 4 .* prop.rc .* r);
+
+ %-- Set up F function for minimization -------------------------------%
+ F = @(r,ii) A1(jj,ii) .* (A2(jj,ii) .* log(C0(jj) ./ m(ii) - ...
+ (sp(jj).omega1 .* r) .^2 )+...
+ A3(jj,ii) .* atanh(sqrt(m(ii) ./ C0(jj)) .* sp(jj).omega1 .* r) + ...
+ A4(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
diff --git a/tfer_diff.m b/tfer_diff.m
new file mode 100644
index 0000000..f009f9e
--- /dev/null
+++ b/tfer_diff.m
@@ -0,0 +1,39 @@
+
+% TFER_DIFF General function for diffusion, only requiring a different G0.
+% Author: Timothy Sipkens, 2018-12-27
+%
+% Inputs:
+% sp Structure defining various setpoint parameters
+% (e.g. m_star, V). Use 'get_setpoint' method to generate
+% this structure.
+% m Particle mass
+% d Particle mobility diameter
+% z Integer charge state
+% prop Device properties (e.g. classifier length)
+%
+% Outputs:
+% Lambda Transfer function
+% G0 Function mapping final to initial radial position
+%=========================================================================%
+
+function [Lambda] = tfer_diff(G0, D, prop)
+
+sig = sqrt(2 .* prop.L .* D ./ prop.v_bar); % diffusive spreading parameter
+
+%-- Evaluate relevant functions ------------------------------------------%
+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 -------------------------%
+% real(...) is added for Case 2C.
+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
+
+end
diff --git a/tfer_ehara.m b/tfer_ehara.m
index 87a8dc2..e74de9e 100644
--- a/tfer_ehara.m
+++ b/tfer_ehara.m
@@ -15,20 +15,23 @@
% Lambda Transfer function
%=========================================================================%
-function [Lambda] = tfer_ehara(sp,m,d,z,prop)
+function [Lambda] = tfer_ehara(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 transfer function -------------------------------------------%
-rho_s = (rs-prop.rc)/prop.del;
-Lambda = ((1-rho_s)+(1+rho_s).*exp(-lam))./2.*and(1m_min).*(m-m_min)./m_del+...
- (m>sp.m_star).*(m m_min) .* (m - m_min) ./ m_del + ...
+ (m > [sp.m_star]') .* (m < [sp.m_max]') .* (([sp.m_star]' - m) ./ m_del + 1);
end
diff --git a/tridiag.m b/tridiag.m
index 2b5df59..7eb2e91 100644
--- a/tridiag.m
+++ b/tridiag.m
@@ -18,7 +18,7 @@
% a(1) is ignored.
%=========================================================================%
-function x = tridiag(a,b,c,y)
+function x = tridiag(a, b, c, y)
%-- Check that the input arrays have acceptable sizes --------------------%
N = length(y);
@@ -38,16 +38,16 @@
y(1) = y(1)/b(1);
for ii = 2:N
- c(ii-1) = c(ii-1)/beta;
- beta = b(ii) - a(ii)*c(ii-1);
- y(ii) = (y(ii) - a(ii)*y(ii-1))/beta;
+ c(ii-1) = c(ii-1) / beta;
+ beta = b(ii) - a(ii) * c(ii-1);
+ y(ii) = (y(ii) - a(ii) * y(ii-1)) / beta;
end
%-- Phase 2: Back-substitution -------------------------------------------%
x(N) = y(N);
for ii = (N-1):-1:1
- x(ii) = y(ii) - c(ii)*x(ii+1);
+ x(ii) = y(ii) - c(ii) * x(ii+1);
end
end