From 9720b8b308544fd3f73bb0b9b8713e55d7444462 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Mon, 15 Jul 2019 13:11:14 -0700 Subject: [PATCH 01/22] Update get_mstar.m --- +tfer_PMA/get_mstar.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m index 5cf143e..aac123b 100644 --- a/+tfer_PMA/get_mstar.m +++ b/+tfer_PMA/get_mstar.m @@ -5,6 +5,8 @@ function [m_star] = get_mstar(prop,V,omega1,omega2) +e = 1.60218e-19; % electron charge [C] + omega_hat = omega2./omega1; alpha = omega1.*(prop.r_hat.^2-omega_hat)./(prop.r_hat.^2-1); From a779c2b3598b62f9cdadb8fd9a3e3e9a286eda3e Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Mon, 15 Jul 2019 13:18:19 -0700 Subject: [PATCH 02/22] Update get_mstar.m --- +tfer_PMA/get_mstar.m | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m index aac123b..787b51d 100644 --- a/+tfer_PMA/get_mstar.m +++ b/+tfer_PMA/get_mstar.m @@ -1,13 +1,20 @@ -% GET_MSTAR Get mass setpoint from provided speeds and voltages +% GET_MSTAR Get mass setpoint for a CPMA from provided speeds and voltages % Author: Timothy Sipkens, 2019-15-07 %=========================================================================% -function [m_star] = get_mstar(prop,V,omega1,omega2) +function [m_star,prop] = get_mstar(prop,V,omega1,omega2) +%-------------------------------------------------------------------------% +% Inputs: +% prop Struct containing physical dimensions of CPMA +% V Operating voltage of the CPMA [V] +% omega1 Rotational speed of inner electrode [rad/s] +% omega2 Rotation speed of outer electrode [rad/s] e = 1.60218e-19; % electron charge [C] omega_hat = omega2./omega1; +prop.omega_hat = omega_hat; % update prop to reflect setpoint alpha = omega1.*(prop.r_hat.^2-omega_hat)./(prop.r_hat.^2-1); beta = omega1.*prop.r1.^2.*(omega_hat-1)./(prop.r_hat^2-1); From a7ab0436b17b6e437d4071cd634257e6bb804f2f Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Mon, 15 Jul 2019 13:25:48 -0700 Subject: [PATCH 03/22] Update get_mstar.m --- +tfer_PMA/get_mstar.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m index 787b51d..584cf63 100644 --- a/+tfer_PMA/get_mstar.m +++ b/+tfer_PMA/get_mstar.m @@ -10,6 +10,12 @@ % V Operating voltage of the CPMA [V] % omega1 Rotational speed of inner electrode [rad/s] % omega2 Rotation speed of outer electrode [rad/s] +% +% Ouputs: +% m_star Mass setpoint [kg] +% prop Updated struct containing physical dimensions of CPMA +% omega_hat is updated +%-------------------------------------------------------------------------% e = 1.60218e-19; % electron charge [C] From a36ee5395b6fc4f9a4740f2957c6304efe0d64f5 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Mon, 15 Jul 2019 13:30:10 -0700 Subject: [PATCH 04/22] Update get_mstar.m --- +tfer_PMA/get_mstar.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m index 584cf63..fe9bc33 100644 --- a/+tfer_PMA/get_mstar.m +++ b/+tfer_PMA/get_mstar.m @@ -13,14 +13,14 @@ % % Ouputs: % m_star Mass setpoint [kg] -% prop Updated struct containing physical dimensions of CPMA +% prop Updated struct containing physical dimensions of CPMA, where % omega_hat is updated %-------------------------------------------------------------------------% e = 1.60218e-19; % electron charge [C] omega_hat = omega2./omega1; -prop.omega_hat = omega_hat; % update prop to reflect setpoint +prop.omega_hat = omega_hat; % update prop to reflect setpoints alpha = omega1.*(prop.r_hat.^2-omega_hat)./(prop.r_hat.^2-1); beta = omega1.*prop.r1.^2.*(omega_hat-1)./(prop.r_hat^2-1); From b4d3b3a3642a917a214af4ea08c5596aa72b51de Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Tue, 16 Jul 2019 21:26:22 -0700 Subject: [PATCH 05/22] Update README.md --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index ca90727..f1d0c0e 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,6 @@ particle tracking methods and using a finite difference method. Information on each file is given as header information in each file, and only a brief overview is provided here. ----------------------------------------------------------------------- ### Code description and components @@ -46,6 +45,7 @@ Note that in these functions, there is a reference to the script (`get_setpoint.m`). This script parses the inputs *d* and *z* and then evaluates the setpoint and related properties, including C0, alpha, and beta. + #### Demonstration script (`main.m`) This script is included to demonstrate evaluation of the transfer function @@ -54,6 +54,7 @@ over multiple cases. Figure 2 that is produced by this procedure will resemble those given in the associated work + #### Remaining functions The remaining functions help in transfer function evaluation, with the @@ -63,10 +64,9 @@ details provided in each file. #### License -This code is distributed under an MIT license -(see corresponding LICENSE file). +This software is licensed under an MIT license (see the corresponding file +for details). ----------------------------------------------------------------------- #### Contact From 98c1dfa9ac27fe90489e28c5e5e91880dc623583 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Wed, 17 Jul 2019 16:06:51 -0700 Subject: [PATCH 06/22] Minor updates to commenting --- +tfer_PMA/get_mstar.m | 2 +- +tfer_PMA/prop_CPMA.m | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m index fe9bc33..0f84ec0 100644 --- a/+tfer_PMA/get_mstar.m +++ b/+tfer_PMA/get_mstar.m @@ -1,5 +1,5 @@ -% GET_MSTAR Get mass setpoint for a CPMA from provided speeds and voltages +% GET_MSTAR Get mass setpoint for a CPMA from provided speeds and voltages. % Author: Timothy Sipkens, 2019-15-07 %=========================================================================% diff --git a/+tfer_PMA/prop_CPMA.m b/+tfer_PMA/prop_CPMA.m index 2fc2634..a473fbb 100644 --- a/+tfer_PMA/prop_CPMA.m +++ b/+tfer_PMA/prop_CPMA.m @@ -1,6 +1,6 @@ % PROP_CPMA Generates the prop struct used to summarize CPMA parameters. -% Author: Timothy Sipkens +% Author: Timothy Sipkens, 2019-06-26 %=========================================================================% function [prop] = prop_CPMA(opt) @@ -15,9 +15,9 @@ if ~exist('opt','var') % if properties set is not specified - opt = 'Buckley'; + opt = 'Olfert'; elseif isempty(opt) - opt = 'Buckley'; + opt = 'Olfert'; end if strcmp(opt,'Olfert') From 4dd253d2cffc54348a726cb07af3be5c2bcc1832 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Wed, 17 Jul 2019 16:36:23 -0700 Subject: [PATCH 07/22] Update README.md --- README.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index f1d0c0e..3c53f0c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -## UBC-tfer-PMA +# UBC-tfer-PMA The attached functions and script are intended to reproduce the results of the associated paper (submitted). They evaluate the transfer function of @@ -49,10 +49,8 @@ evaluates the setpoint and related properties, including C0, alpha, and beta. #### Demonstration script (`main.m`) This script is included to demonstrate evaluation of the transfer function -over multiple cases. - -Figure 2 that is produced by this procedure will resemble those given in -the associated work +over multiple cases. Figure 2 that is produced by this procedure will +resemble those given in the associated work #### Remaining functions From 71fff35cc3447b916ad591033ee0527fddc24586 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Wed, 24 Jul 2019 23:59:53 -0700 Subject: [PATCH 08/22] Removed erroneous values from diff. cases The transfer function output for the various diffusive cases were prone to round-off error and infinite values when using evaluating the error function. These cases are now trimmed out and set to zero. --- +tfer_PMA/tfer_A_diff.m | 4 +++- +tfer_PMA/tfer_B_diff.m | 4 +++- +tfer_PMA/tfer_C_diff.m | 4 +++- +tfer_PMA/tfer_D_diff.m | 4 +++- +tfer_PMA/tfer_E_diff.m | 4 +++- +tfer_PMA/tfer_F_diff.m | 4 +++- 6 files changed, 18 insertions(+), 6 deletions(-) diff --git a/+tfer_PMA/tfer_A_diff.m b/+tfer_PMA/tfer_A_diff.m index c392b3a..ddf16ce 100644 --- a/+tfer_PMA/tfer_A_diff.m +++ b/+tfer_PMA/tfer_A_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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_PMA/tfer_B_diff.m b/+tfer_PMA/tfer_B_diff.m index 4b81396..6d6d025 100644 --- a/+tfer_PMA/tfer_B_diff.m +++ b/+tfer_PMA/tfer_B_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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_PMA/tfer_C_diff.m b/+tfer_PMA/tfer_C_diff.m index f4adb86..9c0e87f 100644 --- a/+tfer_PMA/tfer_C_diff.m +++ b/+tfer_PMA/tfer_C_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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_PMA/tfer_D_diff.m b/+tfer_PMA/tfer_D_diff.m index 0889f79..752da24 100644 --- a/+tfer_PMA/tfer_D_diff.m +++ b/+tfer_PMA/tfer_D_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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_PMA/tfer_E_diff.m b/+tfer_PMA/tfer_E_diff.m index 6e9aa60..04d5397 100644 --- a/+tfer_PMA/tfer_E_diff.m +++ b/+tfer_PMA/tfer_E_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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_PMA/tfer_F_diff.m b/+tfer_PMA/tfer_F_diff.m index afee6e6..6faa458 100644 --- a/+tfer_PMA/tfer_F_diff.m +++ b/+tfer_PMA/tfer_F_diff.m @@ -48,6 +48,8 @@ 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 = max(-1/(4*prop.del).*(K22-K12-K21+K11),0); +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 From 7d384f6761cb11e53cca23e740168e1b459b369e Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 1 Aug 2019 17:04:09 -0700 Subject: [PATCH 09/22] Finite diff. bug fix, +charge script Fixed bug where finite difference fails to evaluate the transfer function for the higher charging states and added a script to evaluate the transfer function over a range of integer charge states. --- +tfer_PMA/tfer_FD.m | 4 +-- main.m | 15 +++++---- main_charge.m | 80 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 90 insertions(+), 9 deletions(-) create mode 100644 main_charge.m diff --git a/+tfer_PMA/tfer_FD.m b/+tfer_PMA/tfer_FD.m index d507723..8f05ada 100644 --- a/+tfer_PMA/tfer_FD.m +++ b/+tfer_PMA/tfer_FD.m @@ -56,8 +56,8 @@ %-- 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>(m_star+2.*sp.m_max),... - m<(m_star-2.*sp.m_max)); + cond0 = or(m>(z.*m_star+2.*sp.m_max),... + m<(z.*m_star-2.*sp.m_max)); % NOTE: conditions limits consideration of those regions where particles % do not escape, speeding computation. ind = ind(~cond0); diff --git a/main.m b/main.m index 7172fee..781f3f2 100644 --- a/main.m +++ b/main.m @@ -10,7 +10,7 @@ Rm = 10; % equivalent resolution of transfer functions (Reavell et al.) -m_star = 1e-18; % mass in kg (1 fg = 1e-18 kg) +m_star = 0.01e-18; % mass in kg (1 fg = 1e-18 kg) m = linspace(0.8,1.2,601).*m_star; % vector of mass z = 1; % integer charge state @@ -148,18 +148,19 @@ m_plot = m./m_star; figure(2); -plot(m_plot,tfer_A); -hold on; +% plot(m_plot,tfer_A); +% hold on; % plot(m_plot,tfer_A_Ehara); -% plot(m_plot,tfer_A_diff); +plot(m_plot,tfer_A_diff); +hold on; % plot(m_plot,tfer_A_pb); -plot(m_plot,tfer_B); +% plot(m_plot,tfer_B); plot(m_plot,tfer_B_diff); % plot(m_plot,tfer_B_pb); % plot(m_plot,tfer_C); -% plot(m_plot,tfer_C_diff); +plot(m_plot,tfer_C_diff); % plot(m_plot,tfer_D); -% plot(m_plot,tfer_D_diff); +plot(m_plot,tfer_D_diff); % plot(m_plot,tfer_E,'r'); % plot(m_plot,tfer_E_diff,'r'); % plot(m_plot,tfer_E_pb); diff --git a/main_charge.m b/main_charge.m new file mode 100644 index 0000000..7faa6a4 --- /dev/null +++ b/main_charge.m @@ -0,0 +1,80 @@ + +% MAIN_CHARGE Script to investigate the effect of multiple chargining. +% Author: Timothy Sipkens, 2019-06-25 +%=========================================================================% + + +%-- Initialize script ----------------------------------------------------% +clear; +close all; + +Rm = 10; % equivalent resolution of transfer functions (Reavell et al.) + +m_star = 0.01e-18; % mass in kg (1 fg = 1e-18 kg) +m = linspace(0.5,4.5,801).*m_star; % vector of mass + +z_max = 4; +for zz=1:z_max + z = zz; % integer charge state + disp(['Processing ',num2str(zz),' of ',num2str(z_max),'...']); + + rho_eff = 900; % effective density + d = (6.*m./(rho_eff.*pi)).^(1/3); + % specify mobility diameter vector with constant effective density + + prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA + % prop.omega_hat = 1; % NOTE: Uncomment for APM condition + + + %=========================================================================% + %-- Finite difference solution -------------------------------------------% + tic; + tfer_FD(:,zz) = tfer_PMA.tfer_FD(m_star,... + m,d,z,prop,'Rm',Rm); + t(1) = toc; + + + %=========================================================================% + %-- Transfer functions for different cases -------------------------------% + %-- Setup for centriputal force ------------------------------------------% + if ~exist('d','var') + B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + else + B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + end + tau = B.*m; + D = prop.D(B).*z; + sig = sqrt(2.*prop.L.*D./prop.v_bar); + D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. + + + %-- Particle tracking approaches -----------------------------------------% + %-- Plug flow ------------------------------------------------------------% + %-- Diffusive transfer functions -----------------------------------------% + %-- Method B --------------------------------% + tic; + tfer_B_diff(:,zz) = tfer_PMA.tfer_B_diff(m_star,m,d,z,prop,'Rm',Rm); + t(12) = toc; +end + +disp('Complete.'); +disp(' '); + + +%=========================================================================% +%-- Plot different transfer functions with respect to m/m* ---------------% +m_plot = m./m_star; + +figure(2); +plot(m_plot,tfer_B_diff); +hold on; +plot(m_plot,min(tfer_FD,1),'k'); +hold off; + +ylim([0,1.2]); +% xlim(1+[-1.5/Rm,1.5/Rm]); + +xlabel('m/m*') +ylabel('{\Lambda}') + + From 0be54ac7977d989d818fa794197e2e3a045251b3 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Fri, 2 Aug 2019 01:50:36 -0700 Subject: [PATCH 10/22] Updates to replicate Ehara --- +tfer_PMA/get_setpoint.m | 36 +++++++++++++---- +tfer_PMA/prop_CPMA.m | 12 +++++- main_Ehara.asv | 85 ++++++++++++++++++++++++++++++++++++++++ main_Ehara.m | 85 ++++++++++++++++++++++++++++++++++++++++ main_charge.m | 19 ++++----- 5 files changed, 218 insertions(+), 19 deletions(-) create mode 100644 main_Ehara.asv create mode 100644 main_Ehara.m diff --git a/+tfer_PMA/get_setpoint.m b/+tfer_PMA/get_setpoint.m index dd3c78c..0f03f78 100644 --- a/+tfer_PMA/get_setpoint.m +++ b/+tfer_PMA/get_setpoint.m @@ -27,10 +27,14 @@ %-------------------------------------------------------------------------% -%-- Set up mobility calculations and parse inputs ------------------------% -if ~exist('z','var') % if integer charge is not specified, use z = 1 - z = 1; -end +%-- Parse inputs ---------------------------------------------------------% +if ~exist('z','var'); z = []; end +if isempty(z); z = 1; end % if integer charge is not specified, use z = 1 + +if ~exist('m_star','var'); m_star = []; end + + +%-- Set up mobility calculations -----------------------------------------% e = 1.60218e-19; % electron charge [C] q = z.*e; % particle charge @@ -48,6 +52,9 @@ if exist('varargin','var') % parse input name-value pair, if it exists if length(varargin)==2 sp.(varargin{1}) = varargin{2}; + elseif length(varargin)==4 + sp.(varargin{1}) = varargin{2}; + sp.(varargin{3}) = varargin{4}; else sp.Rm = 3; % by default use resolution, with value of 3 end @@ -57,18 +64,29 @@ %-- Proceed depnding on which setpoint parameter is specified ------------% -if isfield(sp,'omega1') % if angular speed of inner electrode is specified +if isempty(m_star) % case if m_star is nto specified (use voltage and speed) %-- Azimuth velocity distribution and voltage ------------------------% sp.alpha = sp.omega1.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); - sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; + m_star = sp.V./(log(1/prop.r_hat)./e.*... + (sp.alpha.*prop.rc+sp.beta./prop.rc).^2); + % q = e, z = 1 for setpoint +elseif isfield(sp,'omega1') % if angular speed of inner electrode is specified + + %-- Azimuth velocity distribution and voltage ------------------------% + sp.alpha = sp.omega1.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); + sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); + + sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; + % q = e, z = 1 for setpoint elseif isfield(sp,'V') % if voltage is specified v_theta_rc = sqrt(sp.V.*e./(m_star.*log(1/prop.r_hat))); + % q = e, z = 1 for setpoint A = prop.rc.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1)+... 1./prop.rc.*(prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1)); sp.omega1 = v_theta_rc./A; @@ -82,7 +100,10 @@ %-- Use definition of Rm to derive angular speed at centerline -------% %-- See Reavell et al. (2011) for resolution definition --% n_B = -0.6436; - B_star = tfer_PMA.mp2zp(m_star,1,prop.T,prop.p); % involves invoking mass-mobility relation + B_star = tfer_PMA.mp2zp(m_star,1,prop.T,prop.p); + % involves invoking mass-mobility relation + % z = 1 for the setpoint + sp.m_max = m_star*(1/sp.Rm+1); omega = sqrt(prop.Q/(m_star*B_star*2*pi*prop.rc^2*prop.L*... ((sp.m_max/m_star)^(n_B+1)-(sp.m_max/m_star)^n_B))); @@ -97,7 +118,6 @@ sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; - else error('Invalid setpoint parameter specified.'); diff --git a/+tfer_PMA/prop_CPMA.m b/+tfer_PMA/prop_CPMA.m index a473fbb..5617cc0 100644 --- a/+tfer_PMA/prop_CPMA.m +++ b/+tfer_PMA/prop_CPMA.m @@ -35,13 +35,21 @@ prop.r2 = 0.025; % outer electrode radius [m] prop.r1 = 0.024; % inner electrode radius [m] prop.L = 0.1; % length of APM [m] - RPM = 13350; % rotational speed [rpm] - prop.omega = RPM*2*pi/60; % rotational speed [rad/s] prop.omega_hat = 1; % APM, so rotational speed is the same prop.Q = 1.02e-3/60; % aerosol flowrate [m^3/s] prop.T = 298; % system temperature [K] prop.p = 1; % system pressure [atm] +elseif strcmp(opt,'Ehara') + %-- APM parameters from Ehara et al. -------------% + prop.r2 = 0.103; % outer electrode radius [m] + prop.r1 = 0.1; % inner electrode radius [m] + prop.L = 0.2; % length of APM [m] + prop.omega_hat = 1; % APM, so rotational speed is the same + prop.Q = 1.6667e-5/2; % aerosol flowrate [m^3/s] + prop.T = 298; % system temperature [K] + prop.p = 1; % system pressure [atm] + end diff --git a/main_Ehara.asv b/main_Ehara.asv new file mode 100644 index 0000000..36cdcbd --- /dev/null +++ b/main_Ehara.asv @@ -0,0 +1,85 @@ + +% MAIN Script used in plotting different transfer functions. +% Author: Timothy Sipkens, 2019-06-25 +%=========================================================================% + + +%-- Initialize script ----------------------------------------------------% +clear; +close all; + +V = 10; % voltage to replicate Ehara et al. +omega1 = 3000*0.1047; % angular speed, converted from rpm to rad/s + +e = 1.60218e-19; % electron charge [C] +m = linspace(280,0.45,601)*e; % vector of mass + +z = 1; % integer charge state + +rho_eff = 900; % effective density +d = (6.*m./(rho_eff.*pi)).^(1/3); + % specify mobility diameter vector with constant effective density + +prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA + + +%=========================================================================% +%-- Finite difference solution -------------------------------------------% +tic; +[tfer_FD,~,n] = tfer_PMA.tfer_FD([],... + m,d,1,prop,'V',V,'omega1',omega1); +t(1) = toc; + + +%=========================================================================% +%-- Transfer functions for different cases -------------------------------% +%-- Setup for centriputal force ------------------------------------------% +if ~exist('d','var') + B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); +else + B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +end +tau = B.*m; +D = prop.D(B).*z; +sig = sqrt(2.*prop.L.*D./prop.v_bar); +D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. + + +%-- Particle tracking approaches -----------------------------------------% +%-- Plug flow ------------------------------------------------------------% +%-- Method A ------------------------------% +tic; +[tfer_A,G0_A] = tfer_PMA.tfer_A([],m,d,z,prop,'V',V,'omega1',omega1); +t(2) = toc; + +%-- Method A, Ehara et al. ----------------% +tfer_A_Ehara = tfer_PMA.tfer_A_Ehara([],m,d,z,prop,'V',V,'omega1',omega1); + + + +%-- Parabolic flow -------------------------------------------------------% +%-- Method A ------------------------------% +tic; +[tfer_A_pb,G0_A_pb] = tfer_PMA.tfer_A_pb([],m,d,z,prop,'V',V,'omega1',omega1); +t(8) = toc; + + + +%=========================================================================% +%-- Plot different transfer functions with respect to m/m* ---------------% +m_plot = m./e; + +figure(2); +plot(m_plot,tfer_A); +hold on; +plot(m_plot,tfer_A_Ehara); +plot(m_plot,tfer_A_pb); +% plot(m_plot,min(tfer_FD,1),'k'); +hold off; + +% ylim([0,1.2]); + +xlabel('s') +ylabel('{\Lambda}') + + diff --git a/main_Ehara.m b/main_Ehara.m new file mode 100644 index 0000000..af89e2d --- /dev/null +++ b/main_Ehara.m @@ -0,0 +1,85 @@ + +% MAIN Script used in plotting different transfer functions. +% Author: Timothy Sipkens, 2019-06-25 +%=========================================================================% + + +%-- Initialize script ----------------------------------------------------% +clear; +close all; + +V = 10000; % voltage to replicate Ehara et al. +omega1 = 1000*0.1047; % angular speed, converted from rpm to rad/s + +e = 1.60218e-19; % electron charge [C] +m = linspace(2800,3200,601)*e; % vector of mass + +z = 1; % integer charge state + +rho_eff = 900; % effective density +d = (6.*m./(rho_eff.*pi)).^(1/3); + % specify mobility diameter vector with constant effective density + +prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA + + +%=========================================================================% +%-- Finite difference solution -------------------------------------------% +tic; +[tfer_FD,~,n] = tfer_PMA.tfer_FD([],... + m,d,1,prop,'V',V,'omega1',omega1); +t(1) = toc; + + +%=========================================================================% +%-- Transfer functions for different cases -------------------------------% +%-- Setup for centriputal force ------------------------------------------% +if ~exist('d','var') + B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); +else + B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +end +tau = B.*m; +D = prop.D(B).*z; +sig = sqrt(2.*prop.L.*D./prop.v_bar); +D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. + + +%-- Particle tracking approaches -----------------------------------------% +%-- Plug flow ------------------------------------------------------------% +%-- Method A ------------------------------% +tic; +[tfer_A,G0_A] = tfer_PMA.tfer_A([],m,d,z,prop,'V',V,'omega1',omega1); +t(2) = toc; + +%-- Method A, Ehara et al. ----------------% +tfer_A_Ehara = tfer_PMA.tfer_A_Ehara([],m,d,z,prop,'V',V,'omega1',omega1); + + + +%-- Parabolic flow -------------------------------------------------------% +%-- Method A ------------------------------% +tic; +[tfer_A_pb,G0_A_pb] = tfer_PMA.tfer_A_pb([],m,d,z,prop,'V',V,'omega1',omega1); +t(8) = toc; + + + +%=========================================================================% +%-- Plot different transfer functions with respect to m/m* ---------------% +m_plot = m./e; + +figure(2); +plot(m_plot,tfer_A); +hold on; +plot(m_plot,tfer_A_Ehara); +plot(m_plot,tfer_A_pb); +% plot(m_plot,min(tfer_FD,1),'k'); +hold off; + +% ylim([0,1.2]); + +xlabel('s') +ylabel('{\Lambda}') + + diff --git a/main_charge.m b/main_charge.m index 7faa6a4..ba97ec3 100644 --- a/main_charge.m +++ b/main_charge.m @@ -10,13 +10,14 @@ Rm = 10; % equivalent resolution of transfer functions (Reavell et al.) -m_star = 0.01e-18; % mass in kg (1 fg = 1e-18 kg) -m = linspace(0.5,4.5,801).*m_star; % vector of mass +m_star = 1000e-18; % mass in kg (1 fg = 1e-18 kg) +m = linspace(1e-10,4,801).*m_star; % vector of mass -z_max = 4; -for zz=1:z_max - z = zz; % integer charge state - disp(['Processing ',num2str(zz),' of ',num2str(z_max),'...']); +z_max = 0; +z_vec = 0:z_max; +for zz=1:length(z_vec) + z = z_vec(zz); % integer charge state + disp(['Processing ',num2str(zz),' of ',num2str(length(z_vec)),'...']); rho_eff = 900; % effective density d = (6.*m./(rho_eff.*pi)).^(1/3); @@ -24,12 +25,12 @@ prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA % prop.omega_hat = 1; % NOTE: Uncomment for APM condition - - + + %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; - tfer_FD(:,zz) = tfer_PMA.tfer_FD(m_star,... + [tfer_FD(:,zz),~,n] = tfer_PMA.tfer_FD(m_star,... m,d,z,prop,'Rm',Rm); t(1) = toc; From 1f447c5f987f30c6745cd1a7e0b4d480d054f018 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Fri, 2 Aug 2019 15:14:32 -0700 Subject: [PATCH 11/22] Added gitignore, updates to verification Gitingore is ignoring local copies of +olfert package Uodates to main_Ehara to allow comparison of finite difference --- +tfer_PMA/prop_CPMA.m | 13 ++++++- +tfer_PMA/tfer_FD.m | 29 +++++++++++---- .gitignore | 1 + main_Ehara.asv | 85 ------------------------------------------- main_Ehara.m | 20 +++++----- main_charge.m | 4 +- 6 files changed, 48 insertions(+), 104 deletions(-) create mode 100644 .gitignore delete mode 100644 main_Ehara.asv diff --git a/+tfer_PMA/prop_CPMA.m b/+tfer_PMA/prop_CPMA.m index 5617cc0..05ab72f 100644 --- a/+tfer_PMA/prop_CPMA.m +++ b/+tfer_PMA/prop_CPMA.m @@ -46,10 +46,21 @@ prop.r1 = 0.1; % inner electrode radius [m] prop.L = 0.2; % length of APM [m] prop.omega_hat = 1; % APM, so rotational speed is the same - prop.Q = 1.6667e-5/2; % aerosol flowrate [m^3/s] + prop.Q = 0.5/1000/60; % aerosol flowrate [m^3/s] prop.T = 298; % system temperature [K] prop.p = 1; % system pressure [atm] +elseif strcmp(opt,'Olfert-Collings') + %-- Parameters from Olfert and Collings -------------% + % Nearly identical to the Ehara et al. case + prop.r2 = 0.103; % outer electrode radius [m] + prop.r1 = 0.1; % inner electrode radius [m] + prop.L = 0.2; + prop.omega_hat = 0.945; + prop.Q = 0.5/1000/60; % aerosol flowrate [m^3/s] + prop.T = 295; % system temperature [K] + prop.p = 1; % system pressure [atm] + end diff --git a/+tfer_PMA/tfer_FD.m b/+tfer_PMA/tfer_FD.m index 8f05ada..12724a2 100644 --- a/+tfer_PMA/tfer_FD.m +++ b/+tfer_PMA/tfer_FD.m @@ -75,14 +75,16 @@ ind_mid = 2:(length(r_vec)-1); - %-- Get coefficients ---------------------------------% + + %-- 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); - %-- Righ-hand side of eq. to be solved ---------------% + + %-- 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))+... @@ -92,7 +94,8 @@ (gam-kap(end)+phi(end)).*n(end-1); RHS = @(n) [RHS1(n),RHS2(n),RHS3(n)]; - %-- Form A matrix ------------------------------------% + + %-- Form A matrix -----------------------------------------% b1 = (zet(1)+2*gam+eta(1)); % center c1 = (-gam-kap(1)+phi(1)); % +1 @@ -107,7 +110,8 @@ b = [b1,b2,b3]; c = [c1,c2]; - %-- Initialize variables -----------------------------% + + %-- 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 nargout==3 % initilize variables used for visualizing number concentrations @@ -115,15 +119,26 @@ n_mat(1,:) = n_vec; end - %-- Primary loop for finite difference ---------------% - for jj = 2:nz + + %-- Take first step in implicit scheme, for stability ------% + n_vec = max(tfer_PMA.tridiag([0,a],b,c,RHS(n_vec)),0); + % solve system using Thomas algorithm + + if nargout==3; n_mat(2,:) = n_vec; end + % record particle distribution at this axial position + + + %-- Primary loop for finite difference ---------------------% + for jj = 3:nz n_vec = max(tfer_PMA.tridiag([0,a],b,c,RHS(n_vec)),0); - % solve system using Thoman algorithm + % solve system using Thomas algorithm if nargout==3; n_mat(jj,:) = n_vec; end % record particle distribution at this axial position end + + %-- Format output ------------------------------------------% if nargout==3 % for visualizing number concentrations kk = kk+1; n.n_mat{kk} = max(n_mat,0); diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..13ca4c1 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ ++olfert/ \ No newline at end of file diff --git a/main_Ehara.asv b/main_Ehara.asv deleted file mode 100644 index 36cdcbd..0000000 --- a/main_Ehara.asv +++ /dev/null @@ -1,85 +0,0 @@ - -% MAIN Script used in plotting different transfer functions. -% Author: Timothy Sipkens, 2019-06-25 -%=========================================================================% - - -%-- Initialize script ----------------------------------------------------% -clear; -close all; - -V = 10; % voltage to replicate Ehara et al. -omega1 = 3000*0.1047; % angular speed, converted from rpm to rad/s - -e = 1.60218e-19; % electron charge [C] -m = linspace(280,0.45,601)*e; % vector of mass - -z = 1; % integer charge state - -rho_eff = 900; % effective density -d = (6.*m./(rho_eff.*pi)).^(1/3); - % specify mobility diameter vector with constant effective density - -prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA - - -%=========================================================================% -%-- Finite difference solution -------------------------------------------% -tic; -[tfer_FD,~,n] = tfer_PMA.tfer_FD([],... - m,d,1,prop,'V',V,'omega1',omega1); -t(1) = toc; - - -%=========================================================================% -%-- Transfer functions for different cases -------------------------------% -%-- Setup for centriputal force ------------------------------------------% -if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); -else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); -end -tau = B.*m; -D = prop.D(B).*z; -sig = sqrt(2.*prop.L.*D./prop.v_bar); -D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. - - -%-- Particle tracking approaches -----------------------------------------% -%-- Plug flow ------------------------------------------------------------% -%-- Method A ------------------------------% -tic; -[tfer_A,G0_A] = tfer_PMA.tfer_A([],m,d,z,prop,'V',V,'omega1',omega1); -t(2) = toc; - -%-- Method A, Ehara et al. ----------------% -tfer_A_Ehara = tfer_PMA.tfer_A_Ehara([],m,d,z,prop,'V',V,'omega1',omega1); - - - -%-- Parabolic flow -------------------------------------------------------% -%-- Method A ------------------------------% -tic; -[tfer_A_pb,G0_A_pb] = tfer_PMA.tfer_A_pb([],m,d,z,prop,'V',V,'omega1',omega1); -t(8) = toc; - - - -%=========================================================================% -%-- Plot different transfer functions with respect to m/m* ---------------% -m_plot = m./e; - -figure(2); -plot(m_plot,tfer_A); -hold on; -plot(m_plot,tfer_A_Ehara); -plot(m_plot,tfer_A_pb); -% plot(m_plot,min(tfer_FD,1),'k'); -hold off; - -% ylim([0,1.2]); - -xlabel('s') -ylabel('{\Lambda}') - - diff --git a/main_Ehara.m b/main_Ehara.m index af89e2d..e52a4ad 100644 --- a/main_Ehara.m +++ b/main_Ehara.m @@ -8,11 +8,11 @@ clear; close all; -V = 10000; % voltage to replicate Ehara et al. -omega1 = 1000*0.1047; % angular speed, converted from rpm to rad/s +V = 10; % voltage to replicate Ehara et al. +omega1 = 2500*0.1047; % angular speed, converted from rpm to rad/s e = 1.60218e-19; % electron charge [C] -m = linspace(2800,3200,601)*e; % vector of mass +m = linspace(1e-10,1,601)*e; % vector of mass z = 1; % integer charge state @@ -20,7 +20,9 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA +prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA +% prop.omega_hat = 1; % NOTE: Uncomment for APM condition +% prop.D = @(B) zeros(size(B)); %=========================================================================% @@ -70,11 +72,11 @@ m_plot = m./e; figure(2); -plot(m_plot,tfer_A); -hold on; -plot(m_plot,tfer_A_Ehara); -plot(m_plot,tfer_A_pb); -% plot(m_plot,min(tfer_FD,1),'k'); +% plot(m_plot,tfer_A); +% hold on; +% plot(m_plot,tfer_A_Ehara); +% plot(m_plot,tfer_A_pb); +plot(m_plot,min(tfer_FD,1),'k'); hold off; % ylim([0,1.2]); diff --git a/main_charge.m b/main_charge.m index ba97ec3..a155c7d 100644 --- a/main_charge.m +++ b/main_charge.m @@ -10,10 +10,10 @@ Rm = 10; % equivalent resolution of transfer functions (Reavell et al.) -m_star = 1000e-18; % mass in kg (1 fg = 1e-18 kg) +m_star = 100e-18; % mass in kg (1 fg = 1e-18 kg) m = linspace(1e-10,4,801).*m_star; % vector of mass -z_max = 0; +z_max = 1; z_vec = 0:z_max; for zz=1:length(z_vec) z = z_vec(zz); % integer charge state From fe2859f6e22abb8bff72a782b525bb0770a64020 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Wed, 21 Aug 2019 01:25:47 -0700 Subject: [PATCH 12/22] Case code names, other main scripts - Renamed tfer_*.m functions based on theupdate to case codes. - Added main_charge.m to consider multiple charging - Added main_Ehara, main_kuwata, and main_olfert_collings to replicate figures in other studies - Small update to README - Update to get_setpoint.m to consider the case where m_star is not specified - Update prop_CPMA.m to incorporate other studies - Fixed bug with diffusion --- +tfer_PMA/dm2zp.m | 7 +- +tfer_PMA/get_mstar.m | 32 ------ +tfer_PMA/get_setpoint.m | 67 +++++++---- +tfer_PMA/prop_CPMA.m | 12 +- +tfer_PMA/{tfer_B.m => tfer_1C.m} | 4 +- +tfer_PMA/{tfer_B_diff.m => tfer_1C_diff.m} | 6 +- +tfer_PMA/{tfer_B_pb.m => tfer_1C_pb.m} | 4 +- +tfer_PMA/{tfer_A.m => tfer_1S.m} | 4 +- +tfer_PMA/{tfer_A_diff.m => tfer_1S_diff.m} | 6 +- +tfer_PMA/{tfer_A_pb.m => tfer_1S_pb.m} | 4 +- +tfer_PMA/{tfer_D.m => tfer_2C.m} | 6 +- +tfer_PMA/{tfer_D_diff.m => tfer_2C_diff.m} | 6 +- +tfer_PMA/{tfer_C.m => tfer_2S.m} | 4 +- +tfer_PMA/{tfer_C_diff.m => tfer_2S_diff.m} | 6 +- +tfer_PMA/{tfer_A_Ehara.m => tfer_Ehara.m} | 12 +- +tfer_PMA/tfer_FD.m | 113 ++++++++++-------- +tfer_PMA/{tfer_F.m => tfer_GE.m} | 4 +- +tfer_PMA/{tfer_F_diff.m => tfer_GE_diff.m} | 6 +- +tfer_PMA/{tfer_E.m => tfer_W1.m} | 4 +- +tfer_PMA/{tfer_E_diff.m => tfer_W1_diff.m} | 6 +- +tfer_PMA/{tfer_E_pb.m => tfer_W1_pb.m} | 4 +- README.md | 13 ++- main.m | 121 ++++++++++++-------- main_Ehara.m | 40 +++---- main_charge.m | 31 ++--- main_kuwata.m | 67 +++++++++++ main_olfert_collings.m | 79 +++++++++++++ 27 files changed, 429 insertions(+), 239 deletions(-) delete mode 100644 +tfer_PMA/get_mstar.m rename +tfer_PMA/{tfer_B.m => tfer_1C.m} (92%) rename +tfer_PMA/{tfer_B_diff.m => tfer_1C_diff.m} (91%) rename +tfer_PMA/{tfer_B_pb.m => tfer_1C_pb.m} (94%) rename +tfer_PMA/{tfer_A.m => tfer_1S.m} (93%) rename +tfer_PMA/{tfer_A_diff.m => tfer_1S_diff.m} (90%) rename +tfer_PMA/{tfer_A_pb.m => tfer_1S_pb.m} (94%) rename +tfer_PMA/{tfer_D.m => tfer_2C.m} (90%) rename +tfer_PMA/{tfer_D_diff.m => tfer_2C_diff.m} (91%) rename +tfer_PMA/{tfer_C.m => tfer_2S.m} (93%) rename +tfer_PMA/{tfer_C_diff.m => tfer_2S_diff.m} (91%) rename +tfer_PMA/{tfer_A_Ehara.m => tfer_Ehara.m} (80%) rename +tfer_PMA/{tfer_F.m => tfer_GE.m} (94%) rename +tfer_PMA/{tfer_F_diff.m => tfer_GE_diff.m} (91%) rename +tfer_PMA/{tfer_E.m => tfer_W1.m} (93%) rename +tfer_PMA/{tfer_E_diff.m => tfer_W1_diff.m} (91%) rename +tfer_PMA/{tfer_E_pb.m => tfer_W1_pb.m} (94%) create mode 100644 main_kuwata.m create mode 100644 main_olfert_collings.m diff --git a/+tfer_PMA/dm2zp.m b/+tfer_PMA/dm2zp.m index ab48c1d..ba7adb7 100644 --- a/+tfer_PMA/dm2zp.m +++ b/+tfer_PMA/dm2zp.m @@ -57,11 +57,10 @@ end - +%== CC.m =================================================================% +% Function to evaluate Cunningham slip correction factor. +% Author: Timothy Sipkens, 2019-01-02 function Cc = Cc(d,T,p) -% CC.m Function to evaluate Cunningham slip correction factor. -% Author: Timothy Sipkens, 2019-01-02 -% %-------------------------------------------------------------------------% % Inputs: % d Particle mobility diameter diff --git a/+tfer_PMA/get_mstar.m b/+tfer_PMA/get_mstar.m deleted file mode 100644 index 0f84ec0..0000000 --- a/+tfer_PMA/get_mstar.m +++ /dev/null @@ -1,32 +0,0 @@ - -% GET_MSTAR Get mass setpoint for a CPMA from provided speeds and voltages. -% Author: Timothy Sipkens, 2019-15-07 -%=========================================================================% - -function [m_star,prop] = get_mstar(prop,V,omega1,omega2) -%-------------------------------------------------------------------------% -% Inputs: -% prop Struct containing physical dimensions of CPMA -% V Operating voltage of the CPMA [V] -% omega1 Rotational speed of inner electrode [rad/s] -% omega2 Rotation speed of outer electrode [rad/s] -% -% Ouputs: -% m_star Mass setpoint [kg] -% prop Updated struct containing physical dimensions of CPMA, where -% omega_hat is updated -%-------------------------------------------------------------------------% - -e = 1.60218e-19; % electron charge [C] - -omega_hat = omega2./omega1; -prop.omega_hat = omega_hat; % update prop to reflect setpoints - -alpha = omega1.*(prop.r_hat.^2-omega_hat)./(prop.r_hat.^2-1); -beta = omega1.*prop.r1.^2.*(omega_hat-1)./(prop.r_hat^2-1); - -m_star = V./... - (log(1/prop.r_hat)./e.*(alpha.*prop.rc+beta./prop.rc).^2); - -end - diff --git a/+tfer_PMA/get_setpoint.m b/+tfer_PMA/get_setpoint.m index 0f03f78..164a4f8 100644 --- a/+tfer_PMA/get_setpoint.m +++ b/+tfer_PMA/get_setpoint.m @@ -15,7 +15,7 @@ % ('omega1',double) - Angular speed of inner electrode % ('V',double) - Setpoint voltage % -% Smaple outputs: +% Sample outputs: % C0 Summary parameter for the electrostatic force % tau Product of mechanical mobility and particle mass % sp Struct containing mutliple setpoint parameters (V, alpha, etc.) @@ -34,20 +34,6 @@ if ~exist('m_star','var'); m_star = []; end -%-- Set up mobility calculations -----------------------------------------% -e = 1.60218e-19; % electron charge [C] -q = z.*e; % particle charge - -if ~exist('d','var') % evaluate mechanical mobility - warning('Invoking mass-mobility relation to determine Zp.'); - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); -else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); -end -tau = B.*m; -D = prop.D(B).*z; % diffusion as a function of mechanical mobiltiy and charge state - - %-- Parse inputs for setpoint --------------------------------------------% if exist('varargin','var') % parse input name-value pair, if it exists if length(varargin)==2 @@ -63,17 +49,35 @@ end -%-- Proceed depnding on which setpoint parameter is specified ------------% -if isempty(m_star) % case if m_star is nto specified (use voltage and speed) +%-- Set up mobility calculations -----------------------------------------% +e = 1.60218e-19; % electron charge [C] +q = z.*e; % particle charge + +if ~exist('d','var') % evaluate mechanical mobility + warning('Invoking mass-mobility relation to determine Zp.'); + B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); +else + B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +end +tau = B.*m; +D = prop.D(B).*z; % diffusion as a function of mechanical mobiltiy and charge state + + +%=========================================================================% +%-- Proceed depending on which setpoint parameter is specified -----------% +if isempty(m_star) % case if m_star is not specified (use voltage and speed) %-- Azimuth velocity distribution and voltage ------------------------% - sp.alpha = sp.omega1.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); - sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); + sp.alpha = sp.omega.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); + sp.beta = sp.omega.*prop.rc.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); m_star = sp.V./(log(1/prop.r_hat)./e.*... (sp.alpha.*prop.rc+sp.beta./prop.rc).^2); % q = e, z = 1 for setpoint + sp.omega1 = sp.alpha+sp.beta./(prop.r1.^2); + sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + elseif isfield(sp,'omega1') % if angular speed of inner electrode is specified %-- Azimuth velocity distribution and voltage ------------------------% @@ -83,6 +87,21 @@ sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; % q = e, z = 1 for setpoint + sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + sp.omega = sp.alpha+sp.beta./(prop.rc.^2); + +elseif isfield(sp,'omega') % if angular speed at gap center is specified + + %-- Azimuth velocity distribution and voltage ------------------------% + sp.alpha = sp.omega.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); + sp.beta = sp.omega.*prop.rc.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); + + sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; + % q = e, z = 1 for setpoint + + sp.omega1 = sp.alpha+sp.beta./(prop.r1.^2); + sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + elseif isfield(sp,'V') % if voltage is specified v_theta_rc = sqrt(sp.V.*e./(m_star.*log(1/prop.r_hat))); @@ -93,7 +112,8 @@ sp.alpha = sp.omega1.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); - + sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + sp.omega = sp.alpha+sp.beta./(prop.rc.^2); elseif isfield(sp,'Rm') % if resolution is specified @@ -116,14 +136,19 @@ %-- Azimuth velocity distribution and voltage ------------------------% sp.alpha = sp.omega1.*(prop.r_hat.^2-prop.omega_hat)./(prop.r_hat.^2-1); sp.beta = sp.omega1.*prop.r1.^2.*(prop.omega_hat-1)./(prop.r_hat^2-1); + sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + sp.omega = sp.alpha+sp.beta./(prop.rc.^2); sp.V = m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2; else error('Invalid setpoint parameter specified.'); end +sp.m_star = m_star; + % copy center mass to sp + % center mass is for a singly charged particle - +B_star = tfer_PMA.mp2zp(m_star,1,prop.T,prop.p); C0 = sp.V.*q./log(1/prop.r_hat); % calcualte recurring C0 parameter diff --git a/+tfer_PMA/prop_CPMA.m b/+tfer_PMA/prop_CPMA.m index 05ab72f..5f9c41a 100644 --- a/+tfer_PMA/prop_CPMA.m +++ b/+tfer_PMA/prop_CPMA.m @@ -46,7 +46,7 @@ prop.r1 = 0.1; % inner electrode radius [m] prop.L = 0.2; % length of APM [m] prop.omega_hat = 1; % APM, so rotational speed is the same - prop.Q = 0.5/1000/60; % aerosol flowrate [m^3/s] + prop.Q = 0.5/1000/60; % aerosol flowrate [m^3/s], assumed prop.T = 298; % system temperature [K] prop.p = 1; % system pressure [atm] @@ -61,6 +61,16 @@ prop.T = 295; % system temperature [K] prop.p = 1; % system pressure [atm] +elseif strcmp(opt,'Kuwata') + %-- Parameters from Kuwata --------------------------% + prop.r2 = 0.052; % outer electrode radius [m] + prop.r1 = 0.05; % inner electrode radius [m] + prop.L = 0.25; + prop.omega_hat = 1; + prop.Q = 1.67e-5; % aerosol flowrate [m^3/s] + prop.T = 295; % system temperature [K] + prop.p = 1; % system pressure [atm] + end diff --git a/+tfer_PMA/tfer_B.m b/+tfer_PMA/tfer_1C.m similarity index 92% rename from +tfer_PMA/tfer_B.m rename to +tfer_PMA/tfer_1C.m index b9eb9e0..f1bf7b7 100644 --- a/+tfer_PMA/tfer_B.m +++ b/+tfer_PMA/tfer_1C.m @@ -1,9 +1,9 @@ -% TFER_B Evaluates the transfer function for a PMA in Case B. +% TFER_1C Evaluates the transfer function for a PMA in Case B. % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_B(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1C(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_B_diff.m b/+tfer_PMA/tfer_1C_diff.m similarity index 91% rename from +tfer_PMA/tfer_B_diff.m rename to +tfer_PMA/tfer_1C_diff.m index 6d6d025..7311250 100644 --- a/+tfer_PMA/tfer_B_diff.m +++ b/+tfer_PMA/tfer_1C_diff.m @@ -1,9 +1,9 @@ -% TFER_B_DIFF Evaluates the transfer function for a PMA in Case B (w/ diffusion). +% TFER_1C_DIFF Evaluates the transfer function for a PMA in Case B (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_B_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1C_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_B(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_1C(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_B_pb.m b/+tfer_PMA/tfer_1C_pb.m similarity index 94% rename from +tfer_PMA/tfer_B_pb.m rename to +tfer_PMA/tfer_1C_pb.m index 86901eb..e908f60 100644 --- a/+tfer_PMA/tfer_B_pb.m +++ b/+tfer_PMA/tfer_1C_pb.m @@ -1,9 +1,9 @@ -% TFER_B_PB Evaluates the transfer function for a PMA in Case B (w/ parabolic flow). +% TFER_1C_PB Evaluates the transfer function for a PMA in Case B (w/ parabolic flow). % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_B_pb(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1C_pb(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_A.m b/+tfer_PMA/tfer_1S.m similarity index 93% rename from +tfer_PMA/tfer_A.m rename to +tfer_PMA/tfer_1S.m index 5c0d2c7..163a666 100644 --- a/+tfer_PMA/tfer_A.m +++ b/+tfer_PMA/tfer_1S.m @@ -1,9 +1,9 @@ -% TFER_A Evaluates the transfer function for a PMA in Case A. +% TFER_1S Evaluates the transfer function for a PMA in Case A. % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_A(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1S(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_A_diff.m b/+tfer_PMA/tfer_1S_diff.m similarity index 90% rename from +tfer_PMA/tfer_A_diff.m rename to +tfer_PMA/tfer_1S_diff.m index ddf16ce..67bd157 100644 --- a/+tfer_PMA/tfer_A_diff.m +++ b/+tfer_PMA/tfer_1S_diff.m @@ -1,9 +1,9 @@ -% TFER_A_DIFF Evaluates the transfer function for a PMA in Case A (w/ diffusion). +% TFER_1S_DIFF Evaluates the transfer function for a PMA in Case A (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_A_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1S_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_A(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_1S(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_A_pb.m b/+tfer_PMA/tfer_1S_pb.m similarity index 94% rename from +tfer_PMA/tfer_A_pb.m rename to +tfer_PMA/tfer_1S_pb.m index c41e754..7e433fb 100644 --- a/+tfer_PMA/tfer_A_pb.m +++ b/+tfer_PMA/tfer_1S_pb.m @@ -1,9 +1,9 @@ -% TFER_A_PB Evaluates the transfer function for a PMA in Case A (w/ parabolic flow). +% TFER_1S_PB Evaluates the transfer function for a PMA in Case A (w/ parabolic flow). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_A_pb(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_1S_pb(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_D.m b/+tfer_PMA/tfer_2C.m similarity index 90% rename from +tfer_PMA/tfer_D.m rename to +tfer_PMA/tfer_2C.m index ae192e0..be152d7 100644 --- a/+tfer_PMA/tfer_D.m +++ b/+tfer_PMA/tfer_2C.m @@ -1,9 +1,9 @@ -% TFER_D Evaluates the transfer function for a PMA in Case D. +% TFER_2C Evaluates the transfer function for a PMA in Case D. % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_D(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_2C(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -41,7 +41,7 @@ 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 = real((1/(2*prop.del)).*(rb-ra)); end diff --git a/+tfer_PMA/tfer_D_diff.m b/+tfer_PMA/tfer_2C_diff.m similarity index 91% rename from +tfer_PMA/tfer_D_diff.m rename to +tfer_PMA/tfer_2C_diff.m index 752da24..80384a4 100644 --- a/+tfer_PMA/tfer_D_diff.m +++ b/+tfer_PMA/tfer_2C_diff.m @@ -1,9 +1,9 @@ -% TFER_D_DIFF Evaluates the transfer function for a PMA in Case D (w/ diffusion). +% TFER_2C_DIFF Evaluates the transfer function for a PMA in Case D (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_D_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_2C_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_D(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_2C(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_C.m b/+tfer_PMA/tfer_2S.m similarity index 93% rename from +tfer_PMA/tfer_C.m rename to +tfer_PMA/tfer_2S.m index c7f6662..6de251c 100644 --- a/+tfer_PMA/tfer_C.m +++ b/+tfer_PMA/tfer_2S.m @@ -1,9 +1,9 @@ -% TFER_C Evaluates the transfer function for a PMA in Case C. +% TFER_2S Evaluates the transfer function for a PMA in Case C. % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_C(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_2S(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_C_diff.m b/+tfer_PMA/tfer_2S_diff.m similarity index 91% rename from +tfer_PMA/tfer_C_diff.m rename to +tfer_PMA/tfer_2S_diff.m index 9c0e87f..e4fe62e 100644 --- a/+tfer_PMA/tfer_C_diff.m +++ b/+tfer_PMA/tfer_2S_diff.m @@ -1,9 +1,9 @@ -% TFER_C_DIFF Evaluates the transfer function for a PMA in Case C (w/ diffusion). +% TFER_2S_DIFF Evaluates the transfer function for a PMA in Case C (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_C_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_2S_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_C(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_2S(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_A_Ehara.m b/+tfer_PMA/tfer_Ehara.m similarity index 80% rename from +tfer_PMA/tfer_A_Ehara.m rename to +tfer_PMA/tfer_Ehara.m index 7103847..ef5942b 100644 --- a/+tfer_PMA/tfer_A_Ehara.m +++ b/+tfer_PMA/tfer_Ehara.m @@ -1,9 +1,9 @@ -% TFER_A_EHARA Evaluates the transfer function for a PMA in Case A as per Ehara et al. (1996). +% TFER_EHARA Evaluates the transfer function for a PMA in Case A as per Ehara et al. (1996). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda] = tfer_A_Ehara(m_star,m,d,z,prop,varargin) +function [Lambda] = tfer_Ehara(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -26,9 +26,13 @@ %-- 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 + 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 + rs = real((sqrt(C0./m)+... + sqrt(C0./m-4*sp.alpha*sp.beta))./(2*sp.alpha)); + % equiblirium radius for a given mass end diff --git a/+tfer_PMA/tfer_FD.m b/+tfer_PMA/tfer_FD.m index 12724a2..0eb925c 100644 --- a/+tfer_PMA/tfer_FD.m +++ b/+tfer_PMA/tfer_FD.m @@ -3,7 +3,7 @@ % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,prop,n] = tfer_FD(m_star,m,d,z,prop,varargin) +function [Lambda,sp,n] = tfer_FD(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -44,13 +44,16 @@ del = (prop.r2-prop.r1)/2; rc = (prop.r2+prop.r1)/2; +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); % axial velocity distribution (parabolic) -% v_z = v_bar.*ones(size(r_vec)); % axial velocity distribution (plug) +% v_z = prop.v_bar.*ones(size(r_vec)); % axial velocity distribution (plug) %-- Speed computation using resolution to limit computation --------------% @@ -73,43 +76,6 @@ 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)); - ind_mid = 2:(length(r_vec)-1); - - - %-- 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); - - - %-- 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)]; - - - %-- Form A 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 - - 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]; - %-- Initialize variables -----------------------------------% n_vec = ones(size(r_vec)); % number concentration at current axial position @@ -120,24 +86,31 @@ end - %-- Take first step in implicit scheme, for stability ------% - n_vec = max(tfer_PMA.tridiag([0,a],b,c,RHS(n_vec)),0); - % solve system using Thomas algorithm - - if nargout==3; n_mat(2,:) = n_vec; end - % record particle distribution at this axial position + %-- 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 = 3:nz - n_vec = max(tfer_PMA.tridiag([0,a],b,c,RHS(n_vec)),0); + for jj = 2:nz + n_vec = tfer_PMA.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==3; n_mat(jj,:) = n_vec; end % record particle distribution at this axial position end - %-- Format output ------------------------------------------% if nargout==3 % for visualizing number concentrations kk = kk+1; @@ -153,3 +126,43 @@ end end + + +%== 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) + +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)]; + + +%-- 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 + +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_PMA/tfer_F.m b/+tfer_PMA/tfer_GE.m similarity index 94% rename from +tfer_PMA/tfer_F.m rename to +tfer_PMA/tfer_GE.m index 1900ac0..a30c2e6 100644 --- a/+tfer_PMA/tfer_F.m +++ b/+tfer_PMA/tfer_GE.m @@ -1,9 +1,9 @@ -% TFER_F Evaluates the transfer function for a PMA in Case F. +% TFER_GE Evaluates the transfer function for a PMA in Case F. % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_F(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_GE(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_F_diff.m b/+tfer_PMA/tfer_GE_diff.m similarity index 91% rename from +tfer_PMA/tfer_F_diff.m rename to +tfer_PMA/tfer_GE_diff.m index 6faa458..0d8a6e0 100644 --- a/+tfer_PMA/tfer_F_diff.m +++ b/+tfer_PMA/tfer_GE_diff.m @@ -1,9 +1,9 @@ -% TFER_F_DIFF Evaluates the transfer function for a PMA in Case F (w/ diffusion). +% TFER_GE_DIFF Evaluates the transfer function for a PMA in Case F (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_F_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_GE_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_F(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_GE(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_E.m b/+tfer_PMA/tfer_W1.m similarity index 93% rename from +tfer_PMA/tfer_E.m rename to +tfer_PMA/tfer_W1.m index 4e2aafc..8f806e8 100644 --- a/+tfer_PMA/tfer_E.m +++ b/+tfer_PMA/tfer_W1.m @@ -1,9 +1,9 @@ -% TFER_E Evaluates the transfer function for a PMA in Case E. +% TFER_W1 Evaluates the transfer function for a PMA in Case E. % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_E(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_W1(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass diff --git a/+tfer_PMA/tfer_E_diff.m b/+tfer_PMA/tfer_W1_diff.m similarity index 91% rename from +tfer_PMA/tfer_E_diff.m rename to +tfer_PMA/tfer_W1_diff.m index 04d5397..3c10907 100644 --- a/+tfer_PMA/tfer_E_diff.m +++ b/+tfer_PMA/tfer_W1_diff.m @@ -1,9 +1,9 @@ -% TFER_E_DIFF Evaluates the transfer function for a PMA in Case E (w/ diffusion). +% TFER_W1_DIFF Evaluates the transfer function for a PMA in Case E (w/ diffusion). % Author: Timothy Sipkens, 2018-12-27 %=========================================================================% -function [Lambda,G0] = tfer_E_diff(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_W1_diff(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Inputs: % m_star Setpoint particle mass @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_E(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_PMA.tfer_W1(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_E_pb.m b/+tfer_PMA/tfer_W1_pb.m similarity index 94% rename from +tfer_PMA/tfer_E_pb.m rename to +tfer_PMA/tfer_W1_pb.m index 8e48520..30a05c7 100644 --- a/+tfer_PMA/tfer_E_pb.m +++ b/+tfer_PMA/tfer_W1_pb.m @@ -1,9 +1,9 @@ -% TFER_E_PB Evaluates the transfer function for a PMA in Case E (w/ parabolic flow). +% TFER_W1_PB Evaluates the transfer function for a PMA in Case E (w/ parabolic flow). % Author: Timothy Sipkens, 2019-03-21 %=========================================================================% -function [Lambda,G0] = tfer_E_pb(m_star,m,d,z,prop,varargin) +function [Lambda,G0] = tfer_W1_pb(m_star,m,d,z,prop,varargin) % %-------------------------------------------------------------------------% % Inputs: diff --git a/README.md b/README.md index 3c53f0c..b833a85 100644 --- a/README.md +++ b/README.md @@ -41,6 +41,13 @@ The functions share common inputs: 6. *varargin* (optional) - name-value pairs to specify either the equivalent resolution, inner electrode angular speed, or voltage. +The functions also often share common outputs: + +1. *Lambda* - the transfer function, and + +2. *G0* - the mapping function, transforming a finial radius to the +corresponding position of the particle at the inlet. + Note that in these functions, there is a reference to the script (`get_setpoint.m`). This script parses the inputs *d* and *z* and then evaluates the setpoint and related properties, including C0, alpha, and beta. @@ -50,7 +57,9 @@ evaluates the setpoint and related properties, including C0, alpha, and beta. This script is included to demonstrate evaluation of the transfer function over multiple cases. Figure 2 that is produced by this procedure will -resemble those given in the associated work +resemble those given in the associated work. Other scripts, `main_*.m` +are intended to replicate figures in other works and to consider multiple +charging. #### Remaining functions @@ -70,4 +79,4 @@ for details). This program was written by Timothy A. Sipkens ([tsipkens@mail.ubc.ca](mailto:tsipkens@mail.ubc.ca)) while at the -University of British Columbia +University of British Columbia. diff --git a/main.m b/main.m index 781f3f2..5398874 100644 --- a/main.m +++ b/main.m @@ -20,13 +20,13 @@ % specify mobility diameter vector with constant effective density prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA -% prop.omega_hat = 1; % NOTE: Uncomment for APM condition +prop.omega_hat = 1; % NOTE: Uncomment for APM condition %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; -[tfer_FD,~,n] = tfer_PMA.tfer_FD(m_star,... +[tfer_FD,sp,n] = tfer_PMA.tfer_FD(m_star,... m,d,1,prop,'Rm',Rm); t(1) = toc; @@ -40,99 +40,99 @@ B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); end tau = B.*m; -D = prop.D(B).*z; +D = prop.D(B); sig = sqrt(2.*prop.L.*D./prop.v_bar); D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. %-- Particle tracking approaches -----------------------------------------% %-- Plug flow ------------------------------------------------------------% -%-- Method A ------------------------------% +%-- Method 1S ------------------------------% tic; -[tfer_A,G0_A] = tfer_PMA.tfer_A(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1S,G0_1S] = tfer_PMA.tfer_1S(m_star,m,d,z,prop,'Rm',Rm); t(2) = toc; -%-- Method A, Ehara et al. ----------------% -tfer_A_Ehara = tfer_PMA.tfer_A_Ehara(m_star,m,d,z,prop,'Rm',Rm); +%-- Method 1S, Ehara et al. ----------------% +tfer_Ehara = tfer_PMA.tfer_Ehara(m_star,m,d,z,prop,'Rm',Rm); -%-- Method B ------------------------------% +%-- Method 1C ------------------------------% tic; -[tfer_B,G0_B] = tfer_PMA.tfer_B(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1C,G0_1C] = tfer_PMA.tfer_1C(m_star,m,d,z,prop,'Rm',Rm); t(3) = toc; -%-- Method C ------------------------------% +%-- Method 2S ------------------------------% tic; -[tfer_C,G0_C] = tfer_PMA.tfer_C(m_star,m,d,z,prop,'Rm',Rm); +[tfer_2S,G0_2S] = tfer_PMA.tfer_2S(m_star,m,d,z,prop,'Rm',Rm); t(4) = toc; -%-- Method D ------------------------------% +%-- Method 2C ------------------------------% tic; -[tfer_D,G0_D] = tfer_PMA.tfer_D(m_star,m,d,z,prop,'Rm',Rm); +[tfer_2C,G0_2C] = tfer_PMA.tfer_2C(m_star,m,d,z,prop,'Rm',Rm); t(5) = toc; -%-- Method E ------------------------------% +%-- Method W1 ------------------------------% if prop.omega_hat==1 tic; - [tfer_E,G0_E] = tfer_PMA.tfer_E(m_star,m,d,z,prop,'Rm',Rm); + [tfer_W1,G0_W1] = tfer_PMA.tfer_W1(m_star,m,d,z,prop,'Rm',Rm); t(6) = toc; end -%-- Method F ------------------------------% +%-- Method GE ------------------------------% tic; -[tfer_F,G0_F] = tfer_PMA.tfer_F(m_star,m,d,z,prop,'Rm',Rm); +[tfer_GE,G0_GE] = tfer_PMA.tfer_GE(m_star,m,d,z,prop,'Rm',Rm); t(7) = toc; %-- Parabolic flow -------------------------------------------------------% -%-- Method A ------------------------------% +%-- Method 1S ------------------------------% tic; -[tfer_A_pb,G0_A_pb] = tfer_PMA.tfer_A_pb(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1S_pb,G0_1S_pb] = tfer_PMA.tfer_1S_pb(m_star,m,d,z,prop,'Rm',Rm); t(8) = toc; -%-- Method B ------------------------------% +%-- Method 1C ------------------------------% tic; -[tfer_B_pb,G0_B_pb] = tfer_PMA.tfer_B_pb(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1C_pb,G0_1C_pb] = tfer_PMA.tfer_1C_pb(m_star,m,d,z,prop,'Rm',Rm); t(9) = toc; -%-- Method E ------------------------------% +%-- Method W1 ------------------------------% if prop.omega_hat==1 tic; - [tfer_E_pb,G0_E_pb] = tfer_PMA.tfer_E_pb(m_star,m,d,z,prop,'Rm',Rm); + [tfer_W1_pb,G0_W1_pb] = tfer_PMA.tfer_W1_pb(m_star,m,d,z,prop,'Rm',Rm); t(10) = toc; end %-- Diffusive transfer functions -----------------------------------------% -%-- Method A --------------------------------% +%-- Method 1S --------------------------------% tic; -tfer_A_diff = tfer_PMA.tfer_A_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_1S_diff = tfer_PMA.tfer_1S_diff(m_star,m,d,z,prop,'Rm',Rm); t(11) = toc; -%-- Method B --------------------------------% +%-- Method 1C --------------------------------% tic; -tfer_B_diff = tfer_PMA.tfer_B_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_1C_diff = tfer_PMA.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); t(12) = toc; -%-- Method C --------------------------------% +%-- Method 2S --------------------------------% tic; -tfer_C_diff = tfer_PMA.tfer_C_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_2S_diff = tfer_PMA.tfer_2S_diff(m_star,m,d,z,prop,'Rm',Rm); t(13) = toc; -%-- Method D --------------------------------% +%-- Method 2C --------------------------------% tic; -tfer_D_diff = tfer_PMA.tfer_D_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_2C_diff = tfer_PMA.tfer_2C_diff(m_star,m,d,z,prop,'Rm',Rm); t(14) = toc; -%-- Method E --------------------------------% +%-- Method W1 --------------------------------% if prop.omega_hat==1 tic; - tfer_E_diff = tfer_PMA.tfer_E_diff(m_star,m,d,z,prop,'Rm',Rm); + tfer_W1_diff = tfer_PMA.tfer_W1_diff(m_star,m,d,z,prop,'Rm',Rm); t(15) = toc; end -%-- Method F --------------------------------% +%-- Method GE --------------------------------% tic; -tfer_F_diff = tfer_PMA.tfer_F_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_GE_diff = tfer_PMA.tfer_GE_diff(m_star,m,d,z,prop,'Rm',Rm); t(16) = toc; @@ -148,23 +148,24 @@ m_plot = m./m_star; figure(2); -% plot(m_plot,tfer_A); +% plot(m_plot,tfer_1S); % hold on; -% plot(m_plot,tfer_A_Ehara); -plot(m_plot,tfer_A_diff); +% plot(m_plot,tfer_Ehara); +% plot(m_plot,tfer_1S_diff); +% hold on; +% plot(m_plot,tfer_1S_pb); +plot(m_plot,tfer_1C); hold on; -% plot(m_plot,tfer_A_pb); -% plot(m_plot,tfer_B); -plot(m_plot,tfer_B_diff); -% plot(m_plot,tfer_B_pb); -% plot(m_plot,tfer_C); -plot(m_plot,tfer_C_diff); -% plot(m_plot,tfer_D); -plot(m_plot,tfer_D_diff); -% plot(m_plot,tfer_E,'r'); -% plot(m_plot,tfer_E_diff,'r'); -% plot(m_plot,tfer_E_pb); -% plot(m_plot,tfer_F); +plot(m_plot,tfer_1C_diff); +plot(m_plot,tfer_1C_pb); +% plot(m_plot,tfer_2S); +% plot(m_plot,tfer_2S_diff); +% plot(m_plot,tfer_2C); +% plot(m_plot,tfer_2C_diff); +% plot(m_plot,tfer_W1,'r'); +% plot(m_plot,tfer_W1_diff,'r'); +% plot(m_plot,tfer_W1_pb); +% plot(m_plot,tfer_GE); % plot(m_plot,tfer_tri); plot(m_plot,min(tfer_FD,1),'k'); hold off; @@ -175,4 +176,24 @@ xlabel('m/m*') ylabel('{\Lambda}') +%{ +%=========================================================================% +%-- Bar plot of error ----------------------------------------------------% +chi_sq = []; +mean_sq_err = []; + +vec = {'1S','1C','2S','2C','W1','GE','1S_pb','1C_pb','W1_pb',... + '1S_diff','1C_diff','2S_diff','2C_diff','W1_diff','GE_diff','tri'}; +for ii=1:length(vec) + if ~and(prop.omega_hat~=1,~isempty(strfind(vec{ii},'W'))) + ind_nz = or(eval(['tfer_',vec{ii},'>1e-5']),tfer_FD>1e-5); + chi_sq{ii} = (eval(['tfer_',vec{ii},'(ind_nz)'])-tfer_FD(ind_nz)).^2; + mean_sq_err(ii) = mean(chi_sq{ii}); + end +end + +figure(3); +semilogy(mean_sq_err,'o-'); +ylim([1e-6,1e-2]); +%} diff --git a/main_Ehara.m b/main_Ehara.m index e52a4ad..83b46f1 100644 --- a/main_Ehara.m +++ b/main_Ehara.m @@ -9,10 +9,11 @@ close all; V = 10; % voltage to replicate Ehara et al. -omega1 = 2500*0.1047; % angular speed, converted from rpm to rad/s +omega = 1000*0.1047; % angular speed, converted from rpm to rad/s e = 1.60218e-19; % electron charge [C] -m = linspace(1e-10,1,601)*e; % vector of mass +m = linspace(1e-10,7,601)*e; % vector of mass +% m = linspace(1e-10,6,601)*e; % vector of mass z = 1; % integer charge state @@ -20,49 +21,42 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA -% prop.omega_hat = 1; % NOTE: Uncomment for APM condition -% prop.D = @(B) zeros(size(B)); - +prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; [tfer_FD,~,n] = tfer_PMA.tfer_FD([],... - m,d,1,prop,'V',V,'omega1',omega1); + m,d,1,prop,'V',V,'omega',omega); t(1) = toc; %=========================================================================% %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% -if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); -else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); -end +B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); tau = B.*m; -D = prop.D(B).*z; +D = prop.D(B); sig = sqrt(2.*prop.L.*D./prop.v_bar); D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. %-- Particle tracking approaches -----------------------------------------% %-- Plug flow ------------------------------------------------------------% -%-- Method A ------------------------------% +%-- Method 1S ------------------------------% tic; -[tfer_A,G0_A] = tfer_PMA.tfer_A([],m,d,z,prop,'V',V,'omega1',omega1); +[tfer_1S,G0_1S] = tfer_PMA.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); t(2) = toc; -%-- Method A, Ehara et al. ----------------% -tfer_A_Ehara = tfer_PMA.tfer_A_Ehara([],m,d,z,prop,'V',V,'omega1',omega1); +%-- Method 1S, Ehara et al. ----------------% +tfer_Ehara = tfer_PMA.tfer_Ehara([],m,d,z,prop,'V',V,'omega',omega); %-- Parabolic flow -------------------------------------------------------% -%-- Method A ------------------------------% +%-- Method 1S ------------------------------% tic; -[tfer_A_pb,G0_A_pb] = tfer_PMA.tfer_A_pb([],m,d,z,prop,'V',V,'omega1',omega1); +[tfer_1S_pb,G0_1S_pb] = tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); t(8) = toc; @@ -72,10 +66,10 @@ m_plot = m./e; figure(2); -% plot(m_plot,tfer_A); -% hold on; -% plot(m_plot,tfer_A_Ehara); -% plot(m_plot,tfer_A_pb); +plot(m_plot,tfer_1S); +hold on; +plot(m_plot,tfer_Ehara); +plot(m_plot,tfer_1S_pb); plot(m_plot,min(tfer_FD,1),'k'); hold off; diff --git a/main_charge.m b/main_charge.m index a155c7d..a4b2239 100644 --- a/main_charge.m +++ b/main_charge.m @@ -10,19 +10,19 @@ Rm = 10; % equivalent resolution of transfer functions (Reavell et al.) -m_star = 100e-18; % mass in kg (1 fg = 1e-18 kg) -m = linspace(1e-10,4,801).*m_star; % vector of mass +m_star = 0.01e-18; % mass in kg (1 fg = 1e-18 kg) +m = linspace(1e-10,5,801).*m_star; % vector of mass -z_max = 1; -z_vec = 0:z_max; +z_max = 4; +z_vec = 1:z_max; for zz=1:length(z_vec) z = z_vec(zz); % integer charge state disp(['Processing ',num2str(zz),' of ',num2str(length(z_vec)),'...']); - + rho_eff = 900; % effective density d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density - + prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA % prop.omega_hat = 1; % NOTE: Uncomment for APM condition @@ -30,11 +30,11 @@ %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; - [tfer_FD(:,zz),~,n] = tfer_PMA.tfer_FD(m_star,... + [tfer_FD(:,zz),sp] = tfer_PMA.tfer_FD(m_star,... m,d,z,prop,'Rm',Rm); t(1) = toc; - - + + %=========================================================================% %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% @@ -44,17 +44,18 @@ B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); end tau = B.*m; - D = prop.D(B).*z; + D = prop.D(B); sig = sqrt(2.*prop.L.*D./prop.v_bar); D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. - - + + %-- Particle tracking approaches -----------------------------------------% %-- Plug flow ------------------------------------------------------------% %-- Diffusive transfer functions -----------------------------------------% - %-- Method B --------------------------------% + %-- Method 1C --------------------------------% tic; - tfer_B_diff(:,zz) = tfer_PMA.tfer_B_diff(m_star,m,d,z,prop,'Rm',Rm); + tfer_1C_diff(:,zz) = ... + tfer_PMA.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); t(12) = toc; end @@ -67,7 +68,7 @@ m_plot = m./m_star; figure(2); -plot(m_plot,tfer_B_diff); +plot(m_plot,tfer_1C_diff); hold on; plot(m_plot,min(tfer_FD,1),'k'); hold off; diff --git a/main_kuwata.m b/main_kuwata.m new file mode 100644 index 0000000..009ef32 --- /dev/null +++ b/main_kuwata.m @@ -0,0 +1,67 @@ + +% MAIN Script used in plotting different transfer functions. +% Author: Timothy Sipkens, 2019-06-25 +%=========================================================================% + + +%-- Initialize script ----------------------------------------------------% +clear; +close all; + +V = 100; % voltage to replicate Ehara et al. +omega = 5000*0.1047; % angular speed, converted from rpm to rad/s + +e = 1.60218e-19; % electron charge [C] +m = linspace(0.4,0.7,601).*1e-18; % vector of mass +% m = linspace(1e-10,6,601)*e; % vector of mass + +z = 1; % integer charge state + +d = 100e-9.*ones(size(m)); + % specify mobility diameter vector with constant effective density + +prop = tfer_PMA.prop_CPMA('Kuwata'); % get properties of the CPMA +prop.D = @(B) 1e-10.*ones(size(B)); + + + +%=========================================================================% +%-- Finite difference solution -------------------------------------------% +[tfer_FD,sp] = tfer_PMA.tfer_FD([],... + m,d,1,prop,'V',V,'omega',omega); + + + +%=========================================================================% +%-- Transfer functions for different cases -------------------------------% +%-- Setup for centriputal force ------------------------------------------% +B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +tau = B.*m; + +%-- Particle tracking approaches -----------------------------------------% +%-- Plug flow ------------------------------------------------------------% +%-- Method 1S ------------------------------% +[tfer_1S] = ... + tfer_PMA.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); +[tfer_1S_pb] = ... + tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + + + +%=========================================================================% +%-- Plot different transfer functions with respect to m/m* ---------------% +m_plot = m; + +figure(2); +plot(m_plot,min(tfer_FD,1)); +hold on; +plot(m_plot,min(tfer_1S,1)); +plot(m_plot,min(tfer_1S_pb,1)); +hold off; + +% ylim([0,1.2]); + +xlabel('s') +ylabel('{\Lambda}') + + diff --git a/main_olfert_collings.m b/main_olfert_collings.m new file mode 100644 index 0000000..a3abab9 --- /dev/null +++ b/main_olfert_collings.m @@ -0,0 +1,79 @@ + +% MAIN Script used in plotting different transfer functions. +% Author: Timothy Sipkens, 2019-06-25 +%=========================================================================% + + +%-- Initialize script ----------------------------------------------------% +clear; +close all; + +V = 10; % voltage to replicate Ehara et al. +omega = 2500*0.1047; % angular speed, converted from rpm to rad/s + +e = 1.60218e-19; % electron charge [C] +m = linspace(1e-10,1,601)*e; % vector of mass +% m = linspace(1e-10,6,601)*e; % vector of mass + +z = 1; % integer charge state + +rho_eff = 900; % effective density +d = (6.*m./(rho_eff.*pi)).^(1/3); + % specify mobility diameter vector with constant effective density + +prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA +prop.D = @(B) 1e-10.*ones(size(B)); +omega_hat = prop.omega_hat; % only valid for CPMA + + + +%=========================================================================% +%-- Finite difference solution -------------------------------------------% +prop.omega_hat = 1; +[tfer_FD_w1,sp] = tfer_PMA.tfer_FD([],... + m,d,1,prop,'V',V,'omega',omega); + +prop.omega_hat = omega_hat; +[tfer_FD,sp_cpma] = tfer_PMA.tfer_FD([],... + m,d,1,prop,'V',V,'omega',omega); + + + +%=========================================================================% +%-- Transfer functions for different cases -------------------------------% +%-- Setup for centriputal force ------------------------------------------% +prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA +B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +tau = B.*m; + +%-- Particle tracking approaches -----------------------------------------% +%-- Plug flow ------------------------------------------------------------% +%-- Method 1S ------------------------------% +prop.omega_hat = 1; % NOTE: Uncomment for APM condition +[tfer_1S_w1] = ... + tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + +prop.omega_hat = omega_hat; +[tfer_1S] = ... + tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + + + +%=========================================================================% +%-- Plot different transfer functions with respect to m/m* ---------------% +m_plot = m./e; + +figure(2); +plot(m_plot,min(tfer_FD,1)); +hold on; +plot(m_plot,min(tfer_FD_w1,1)); +plot(m_plot,min(tfer_1S,1)); +plot(m_plot,min(tfer_1S_w1,1)); +hold off; + +% ylim([0,1.2]); + +xlabel('s') +ylabel('{\Lambda}') + + From 52da9b6a9b3caccbf69e6877ee2fa78601ac3f8c Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Wed, 21 Aug 2019 22:16:07 -0700 Subject: [PATCH 13/22] Formatting and working updates --- +tfer_PMA/G_fun.m | 27 ++++++++++++--------------- +tfer_PMA/tfer_FD.m | 4 +--- main_charge.m | 4 ++-- 3 files changed, 15 insertions(+), 20 deletions(-) diff --git a/+tfer_PMA/G_fun.m b/+tfer_PMA/G_fun.m index 0193705..7bc7d94 100644 --- a/+tfer_PMA/G_fun.m +++ b/+tfer_PMA/G_fun.m @@ -19,7 +19,8 @@ %-------------------------------------------------------------------------% -condit = (alpha^2) < (beta^2./(rs.^4)); % determines where on is relative to rs +condit = (alpha^2) < (beta^2./(rs.^4)); + % whether or not lambda is positive of negative ns = length(rs); % number of equilirbrium radii to consider (one per particle mass considered) nL = length(rL); % number of final radial positions to consider (usually only one) @@ -33,7 +34,7 @@ for jj=1:nL % loop throught all specified particle exit radii for ii=1:ns % loop to consider multiple equilbirium radii (multiple m) - if rL(jj) (rs(ii)-5*offset)) - % if rL = rs, then r0 = rs + % if rL = rs, within offset, then r0 = rs G(jj,ii) = rL(jj); elseif rL(jj)>rs(ii) % if rL > rs, then r0 > rL - if condit(ii) + if condit(ii) % zero is in interval [rL,r2] if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),r2,ii)) % if sign does not change in interval [rL,r2], return r2 G(jj,ii) = r2; - else - % find zero in interval [rL,r2] + else % find zero in interval [rL,r2] G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [rL(jj),r2]); end - else + else % zero is in interval [max(r1,rs),rL] if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),max(r1,rs(ii)+offset),ii)) % if sign does not change in interval [max(r1,rs),rL], return max(r1,rs) G(jj,ii) = max(r1,rs(ii)+offset); - else - % find zero in interval [max(r1,rs),rL] + else % find zero in interval [max(r1,rs),rL] G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [max(r1,rs(ii)+offset),rL(jj)]); end @@ -74,24 +73,22 @@ elseif rL(jj) Date: Thu, 12 Sep 2019 10:46:56 -0700 Subject: [PATCH 14/22] Updated name of prop_PMA --- +tfer_PMA/{prop_CPMA.m => prop_PMA.m} | 0 main.m | 2 +- main_Ehara.m | 2 +- main_charge.m | 2 +- main_kuwata.m | 2 +- main_olfert_collings.m | 2 +- 6 files changed, 5 insertions(+), 5 deletions(-) rename +tfer_PMA/{prop_CPMA.m => prop_PMA.m} (100%) diff --git a/+tfer_PMA/prop_CPMA.m b/+tfer_PMA/prop_PMA.m similarity index 100% rename from +tfer_PMA/prop_CPMA.m rename to +tfer_PMA/prop_PMA.m diff --git a/main.m b/main.m index 5398874..15ab633 100644 --- a/main.m +++ b/main.m @@ -19,7 +19,7 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA +prop = tfer_PMA.prop_PMA('Olfert'); % get properties of the CPMA prop.omega_hat = 1; % NOTE: Uncomment for APM condition diff --git a/main_Ehara.m b/main_Ehara.m index 83b46f1..4dec54b 100644 --- a/main_Ehara.m +++ b/main_Ehara.m @@ -21,7 +21,7 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Ehara'); % get properties of the CPMA +prop = tfer_PMA.prop_PMA('Ehara'); % get properties of the CPMA %=========================================================================% %-- Finite difference solution -------------------------------------------% diff --git a/main_charge.m b/main_charge.m index aceba88..3840786 100644 --- a/main_charge.m +++ b/main_charge.m @@ -23,7 +23,7 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density - prop = tfer_PMA.prop_CPMA('Olfert'); % get properties of the CPMA + prop = tfer_PMA.prop_PMA('Olfert'); % get properties of the CPMA % prop.omega_hat = 1; % NOTE: Uncomment for APM condition diff --git a/main_kuwata.m b/main_kuwata.m index 009ef32..47b08ee 100644 --- a/main_kuwata.m +++ b/main_kuwata.m @@ -20,7 +20,7 @@ d = 100e-9.*ones(size(m)); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Kuwata'); % get properties of the CPMA +prop = tfer_PMA.prop_PMA('Kuwata'); % get properties of the CPMA prop.D = @(B) 1e-10.*ones(size(B)); diff --git a/main_olfert_collings.m b/main_olfert_collings.m index a3abab9..92550fd 100644 --- a/main_olfert_collings.m +++ b/main_olfert_collings.m @@ -21,7 +21,7 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA +prop = tfer_PMA.prop_PMA('Olfert-Collings'); % get properties of the CPMA prop.D = @(B) 1e-10.*ones(size(B)); omega_hat = prop.omega_hat; % only valid for CPMA From 2ff40b774793b968c0f15a384614a4e2a4fb91ba Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 19 Sep 2019 15:18:53 -0700 Subject: [PATCH 15/22] Updates to naming convention --- +tfer_PMA/get_setpoint.m | 21 ++++++++++++++------ +tfer_PMA/mp2zp.m | 8 ++++---- +tfer_PMA/prop_PMA.m | 2 +- +tfer_PMA/tfer_1C.m | 2 +- +tfer_PMA/tfer_1C_diff.m | 6 +++--- +tfer_PMA/tfer_1C_pb.m | 4 ++-- +tfer_PMA/tfer_1S.m | 2 +- +tfer_PMA/tfer_1S_diff.m | 6 +++--- +tfer_PMA/tfer_1S_pb.m | 4 ++-- +tfer_PMA/tfer_2C.m | 2 +- +tfer_PMA/tfer_2C_diff.m | 6 +++--- +tfer_PMA/tfer_2S.m | 2 +- +tfer_PMA/tfer_2S_diff.m | 6 +++--- +tfer_PMA/tfer_Ehara.m | 2 +- +tfer_PMA/tfer_FD.m | 4 ++-- +tfer_PMA/tfer_GE.m | 4 ++-- +tfer_PMA/tfer_GE_diff.m | 6 +++--- +tfer_PMA/tfer_W1.m | 2 +- +tfer_PMA/tfer_W1_diff.m | 6 +++--- +tfer_PMA/tfer_W1_pb.m | 4 ++-- +tfer_PMA/tfer_tri.m | 2 +- main.m | 42 ++++++++++++++++++++-------------------- main_Ehara.m | 12 ++++++------ main_charge.m | 10 +++++----- main_kuwata.m | 10 +++++----- main_olfert_collings.m | 14 +++++++------- 26 files changed, 99 insertions(+), 90 deletions(-) diff --git a/+tfer_PMA/get_setpoint.m b/+tfer_PMA/get_setpoint.m index 164a4f8..80e6044 100644 --- a/+tfer_PMA/get_setpoint.m +++ b/+tfer_PMA/get_setpoint.m @@ -55,9 +55,9 @@ if ~exist('d','var') % evaluate mechanical mobility warning('Invoking mass-mobility relation to determine Zp.'); - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end tau = B.*m; D = prop.D(B).*z; % diffusion as a function of mechanical mobiltiy and charge state @@ -78,6 +78,15 @@ sp.omega1 = sp.alpha+sp.beta./(prop.r1.^2); sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); + %-- Calculate resolution ---------------------------------------------% + n_B = -0.6436; + B_star = tfer_pma.mp2zp(m_star,1,prop.T,prop.p); + t0 = prop.Q/(m_star*B_star*2*pi*prop.L*... + sp.omega^2*prop.rc^2); + m_rat = @(Rm) 1/Rm+1; + fun = @(Rm) (m_rat(Rm))^(n_B+1)-(m_rat(Rm))^n_B; + sp.Rm = fminsearch(@(Rm) (t0-fun(Rm))^2,10); + elseif isfield(sp,'omega1') % if angular speed of inner electrode is specified %-- Azimuth velocity distribution and voltage ------------------------% @@ -89,7 +98,7 @@ sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2); sp.omega = sp.alpha+sp.beta./(prop.rc.^2); - + elseif isfield(sp,'omega') % if angular speed at gap center is specified %-- Azimuth velocity distribution and voltage ------------------------% @@ -120,10 +129,10 @@ %-- Use definition of Rm to derive angular speed at centerline -------% %-- See Reavell et al. (2011) for resolution definition --% n_B = -0.6436; - B_star = tfer_PMA.mp2zp(m_star,1,prop.T,prop.p); + B_star = tfer_pma.mp2zp(m_star,1,prop.T,prop.p); % involves invoking mass-mobility relation % z = 1 for the setpoint - + sp.m_max = m_star*(1/sp.Rm+1); omega = sqrt(prop.Q/(m_star*B_star*2*pi*prop.rc^2*prop.L*... ((sp.m_max/m_star)^(n_B+1)-(sp.m_max/m_star)^n_B))); @@ -148,7 +157,7 @@ % copy center mass to sp % center mass is for a singly charged particle -B_star = tfer_PMA.mp2zp(m_star,1,prop.T,prop.p); +B_star = tfer_pma.mp2zp(m_star,1,prop.T,prop.p); C0 = sp.V.*q./log(1/prop.r_hat); % calcualte recurring C0 parameter diff --git a/+tfer_PMA/mp2zp.m b/+tfer_PMA/mp2zp.m index 49ba9ac..f928c07 100644 --- a/+tfer_PMA/mp2zp.m +++ b/+tfer_PMA/mp2zp.m @@ -24,17 +24,17 @@ %-- Invoke mass-mobility relation ----------------------------------------% -mass_mob_pref = 524; -mass_mob_exp = 3; +mass_mob_pref = 0.0612; %524; +mass_mob_exp = 2.48; %3; d = (m./mass_mob_pref).^(1/mass_mob_exp); % use mass-mobility relationship to get mobility diameter %-- Use mobility diameter to get particle electro and mechanical mobl. ---% if nargin<3 - [Zp,B] = tfer_PMA.dm2zp(d,z); + [Zp,B] = tfer_pma.dm2zp(d,z); else - [Zp,B] = tfer_PMA.dm2zp(d,z,T,P); + [Zp,B] = tfer_pma.dm2zp(d,z,T,P); end end diff --git a/+tfer_PMA/prop_PMA.m b/+tfer_PMA/prop_PMA.m index 5f9c41a..3b6ed29 100644 --- a/+tfer_PMA/prop_PMA.m +++ b/+tfer_PMA/prop_PMA.m @@ -3,7 +3,7 @@ % Author: Timothy Sipkens, 2019-06-26 %=========================================================================% -function [prop] = prop_CPMA(opt) +function [prop] = prop_PMA(opt) %-------------------------------------------------------------------------% % Input: % opt Options string specifying parameter set diff --git a/+tfer_PMA/tfer_1C.m b/+tfer_PMA/tfer_1C.m index f1bf7b7..ec4ab2b 100644 --- a/+tfer_PMA/tfer_1C.m +++ b/+tfer_PMA/tfer_1C.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- 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)); diff --git a/+tfer_PMA/tfer_1C_diff.m b/+tfer_PMA/tfer_1C_diff.m index 7311250..7867aa0 100644 --- a/+tfer_PMA/tfer_1C_diff.m +++ b/+tfer_PMA/tfer_1C_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_1C(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_1C(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_1C_pb.m b/+tfer_PMA/tfer_1C_pb.m index e908f60..233c30b 100644 --- a/+tfer_PMA/tfer_1C_pb.m +++ b/+tfer_PMA/tfer_1C_pb.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- 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)); @@ -47,7 +47,7 @@ %-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) tfer_PMA.G_fun(min_fun,r,rs,prop.r1,prop.r2,sp.alpha,sp.beta); +G0 = @(r) tfer_pma.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))); diff --git a/+tfer_PMA/tfer_1S.m b/+tfer_PMA/tfer_1S.m index 163a666..caaa97e 100644 --- a/+tfer_PMA/tfer_1S.m +++ b/+tfer_PMA/tfer_1S.m @@ -21,7 +21,7 @@ % G0 Function mapping final to initial radial position %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_PMA/tfer_1S_diff.m b/+tfer_PMA/tfer_1S_diff.m index 67bd157..fadcb96 100644 --- a/+tfer_PMA/tfer_1S_diff.m +++ b/+tfer_PMA/tfer_1S_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_1S(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_1S(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_1S_pb.m b/+tfer_PMA/tfer_1S_pb.m index 7e433fb..d0eefdf 100644 --- a/+tfer_PMA/tfer_1S_pb.m +++ b/+tfer_PMA/tfer_1S_pb.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc @@ -46,7 +46,7 @@ %-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) tfer_PMA.G_fun(min_fun,r,rs,prop.r1,prop.r2,sp.alpha,sp.beta); +G0 = @(r) tfer_pma.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))); diff --git a/+tfer_PMA/tfer_2C.m b/+tfer_PMA/tfer_2C.m index be152d7..0b75275 100644 --- a/+tfer_PMA/tfer_2C.m +++ b/+tfer_PMA/tfer_2C.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- 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)); diff --git a/+tfer_PMA/tfer_2C_diff.m b/+tfer_PMA/tfer_2C_diff.m index 80384a4..97120e3 100644 --- a/+tfer_PMA/tfer_2C_diff.m +++ b/+tfer_PMA/tfer_2C_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_2C(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_2C(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_2S.m b/+tfer_PMA/tfer_2S.m index 6de251c..06ecc44 100644 --- a/+tfer_PMA/tfer_2S.m +++ b/+tfer_PMA/tfer_2S.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_PMA/tfer_2S_diff.m b/+tfer_PMA/tfer_2S_diff.m index e4fe62e..9045298 100644 --- a/+tfer_PMA/tfer_2S_diff.m +++ b/+tfer_PMA/tfer_2S_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_2S(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_2S(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_Ehara.m b/+tfer_PMA/tfer_Ehara.m index ef5942b..f880d71 100644 --- a/+tfer_PMA/tfer_Ehara.m +++ b/+tfer_PMA/tfer_Ehara.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_PMA/tfer_FD.m b/+tfer_PMA/tfer_FD.m index 1cb1591..aa077d1 100644 --- a/+tfer_PMA/tfer_FD.m +++ b/+tfer_PMA/tfer_FD.m @@ -27,7 +27,7 @@ % laboratory. %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Discretize the space -------------------------------------------------% @@ -97,7 +97,7 @@ %-- Primary loop for finite difference ---------------------% for jj = 2:nz - n_vec = tfer_PMA.tridiag([0,a],b,c,RHS(n_vec)); + n_vec = tfer_pma.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 diff --git a/+tfer_PMA/tfer_GE.m b/+tfer_PMA/tfer_GE.m index a30c2e6..9d0b3fe 100644 --- a/+tfer_PMA/tfer_GE.m +++ b/+tfer_PMA/tfer_GE.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc @@ -47,7 +47,7 @@ %-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) tfer_PMA.G_fun(min_fun,r,rs,prop.r1,prop.r2,sp.alpha,sp.beta); +G0 = @(r) tfer_pma.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))); diff --git a/+tfer_PMA/tfer_GE_diff.m b/+tfer_PMA/tfer_GE_diff.m index 0d8a6e0..015974b 100644 --- a/+tfer_PMA/tfer_GE_diff.m +++ b/+tfer_PMA/tfer_GE_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_GE(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_GE(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_W1.m b/+tfer_PMA/tfer_W1.m index 8f806e8..b9751a7 100644 --- a/+tfer_PMA/tfer_W1.m +++ b/+tfer_PMA/tfer_W1.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_PMA/tfer_W1_diff.m b/+tfer_PMA/tfer_W1_diff.m index 3c10907..d3741eb 100644 --- a/+tfer_PMA/tfer_W1_diff.m +++ b/+tfer_PMA/tfer_W1_diff.m @@ -24,10 +24,10 @@ %-- Evaluate mechanical mobility for diffusion calc. ---------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); % if mobility is not specified, use mass-mobility relation to estimate else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end D = prop.D(B).*z; @@ -35,7 +35,7 @@ % integer charge state sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter -[~,G0] = tfer_PMA.tfer_W1(m_star,m,d,z,prop,varargin{:}); +[~,G0] = tfer_pma.tfer_W1(m_star,m,d,z,prop,varargin{:}); % get G0 function for this case rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring quantity diff --git a/+tfer_PMA/tfer_W1_pb.m b/+tfer_PMA/tfer_W1_pb.m index 30a05c7..b84ecc6 100644 --- a/+tfer_PMA/tfer_W1_pb.m +++ b/+tfer_PMA/tfer_W1_pb.m @@ -23,7 +23,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint (parses d and z) +tfer_pma.get_setpoint; % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc @@ -47,7 +47,7 @@ %-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) tfer_PMA.G_fun(min_fun,r,rs,prop.r1,prop.r2,sp.alpha,sp.beta); +G0 = @(r) tfer_pma.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))); diff --git a/+tfer_PMA/tfer_tri.m b/+tfer_PMA/tfer_tri.m index b743af3..69ca2da 100644 --- a/+tfer_PMA/tfer_tri.m +++ b/+tfer_PMA/tfer_tri.m @@ -22,7 +22,7 @@ %-------------------------------------------------------------------------% -tfer_PMA.get_setpoint; % get setpoint +tfer_pma.get_setpoint; % get setpoint if ~isfield(sp,'m_max') % if m_max was not specified n_B = -0.6436; diff --git a/main.m b/main.m index 15ab633..445e268 100644 --- a/main.m +++ b/main.m @@ -19,14 +19,14 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_PMA('Olfert'); % get properties of the CPMA +prop = tfer_pma.prop_PMA('Olfert'); % get properties of the CPMA prop.omega_hat = 1; % NOTE: Uncomment for APM condition %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; -[tfer_FD,sp,n] = tfer_PMA.tfer_FD(m_star,... +[tfer_FD,sp,n] = tfer_pma.tfer_FD(m_star,... m,d,1,prop,'Rm',Rm); t(1) = toc; @@ -35,9 +35,9 @@ %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end tau = B.*m; D = prop.D(B); @@ -49,55 +49,55 @@ %-- Plug flow ------------------------------------------------------------% %-- Method 1S ------------------------------% tic; -[tfer_1S,G0_1S] = tfer_PMA.tfer_1S(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1S,G0_1S] = tfer_pma.tfer_1S(m_star,m,d,z,prop,'Rm',Rm); t(2) = toc; %-- Method 1S, Ehara et al. ----------------% -tfer_Ehara = tfer_PMA.tfer_Ehara(m_star,m,d,z,prop,'Rm',Rm); +tfer_Ehara = tfer_pma.tfer_Ehara(m_star,m,d,z,prop,'Rm',Rm); %-- Method 1C ------------------------------% tic; -[tfer_1C,G0_1C] = tfer_PMA.tfer_1C(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1C,G0_1C] = tfer_pma.tfer_1C(m_star,m,d,z,prop,'Rm',Rm); t(3) = toc; %-- Method 2S ------------------------------% tic; -[tfer_2S,G0_2S] = tfer_PMA.tfer_2S(m_star,m,d,z,prop,'Rm',Rm); +[tfer_2S,G0_2S] = tfer_pma.tfer_2S(m_star,m,d,z,prop,'Rm',Rm); t(4) = toc; %-- Method 2C ------------------------------% tic; -[tfer_2C,G0_2C] = tfer_PMA.tfer_2C(m_star,m,d,z,prop,'Rm',Rm); +[tfer_2C,G0_2C] = tfer_pma.tfer_2C(m_star,m,d,z,prop,'Rm',Rm); t(5) = toc; %-- Method W1 ------------------------------% if prop.omega_hat==1 tic; - [tfer_W1,G0_W1] = tfer_PMA.tfer_W1(m_star,m,d,z,prop,'Rm',Rm); + [tfer_W1,G0_W1] = tfer_pma.tfer_W1(m_star,m,d,z,prop,'Rm',Rm); t(6) = toc; end %-- Method GE ------------------------------% tic; -[tfer_GE,G0_GE] = tfer_PMA.tfer_GE(m_star,m,d,z,prop,'Rm',Rm); +[tfer_GE,G0_GE] = tfer_pma.tfer_GE(m_star,m,d,z,prop,'Rm',Rm); t(7) = toc; %-- Parabolic flow -------------------------------------------------------% %-- Method 1S ------------------------------% tic; -[tfer_1S_pb,G0_1S_pb] = tfer_PMA.tfer_1S_pb(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1S_pb,G0_1S_pb] = tfer_pma.tfer_1S_pb(m_star,m,d,z,prop,'Rm',Rm); t(8) = toc; %-- Method 1C ------------------------------% tic; -[tfer_1C_pb,G0_1C_pb] = tfer_PMA.tfer_1C_pb(m_star,m,d,z,prop,'Rm',Rm); +[tfer_1C_pb,G0_1C_pb] = tfer_pma.tfer_1C_pb(m_star,m,d,z,prop,'Rm',Rm); t(9) = toc; %-- Method W1 ------------------------------% if prop.omega_hat==1 tic; - [tfer_W1_pb,G0_W1_pb] = tfer_PMA.tfer_W1_pb(m_star,m,d,z,prop,'Rm',Rm); + [tfer_W1_pb,G0_W1_pb] = tfer_pma.tfer_W1_pb(m_star,m,d,z,prop,'Rm',Rm); t(10) = toc; end @@ -105,40 +105,40 @@ %-- Diffusive transfer functions -----------------------------------------% %-- Method 1S --------------------------------% tic; -tfer_1S_diff = tfer_PMA.tfer_1S_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_1S_diff = tfer_pma.tfer_1S_diff(m_star,m,d,z,prop,'Rm',Rm); t(11) = toc; %-- Method 1C --------------------------------% tic; -tfer_1C_diff = tfer_PMA.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_1C_diff = tfer_pma.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); t(12) = toc; %-- Method 2S --------------------------------% tic; -tfer_2S_diff = tfer_PMA.tfer_2S_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_2S_diff = tfer_pma.tfer_2S_diff(m_star,m,d,z,prop,'Rm',Rm); t(13) = toc; %-- Method 2C --------------------------------% tic; -tfer_2C_diff = tfer_PMA.tfer_2C_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_2C_diff = tfer_pma.tfer_2C_diff(m_star,m,d,z,prop,'Rm',Rm); t(14) = toc; %-- Method W1 --------------------------------% if prop.omega_hat==1 tic; - tfer_W1_diff = tfer_PMA.tfer_W1_diff(m_star,m,d,z,prop,'Rm',Rm); + tfer_W1_diff = tfer_pma.tfer_W1_diff(m_star,m,d,z,prop,'Rm',Rm); t(15) = toc; end %-- Method GE --------------------------------% tic; -tfer_GE_diff = tfer_PMA.tfer_GE_diff(m_star,m,d,z,prop,'Rm',Rm); +tfer_GE_diff = tfer_pma.tfer_GE_diff(m_star,m,d,z,prop,'Rm',Rm); t(16) = toc; %-- Triangle approx. -----------------------% tic; -tfer_tri = tfer_PMA.tfer_tri(m_star,m,d,z,prop,'Rm',Rm); +tfer_tri = tfer_pma.tfer_tri(m_star,m,d,z,prop,'Rm',Rm); t(18) = toc; diff --git a/main_Ehara.m b/main_Ehara.m index 4dec54b..0482bb0 100644 --- a/main_Ehara.m +++ b/main_Ehara.m @@ -21,12 +21,12 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_PMA('Ehara'); % get properties of the CPMA +prop = tfer_pma.prop_PMA('Ehara'); % get properties of the CPMA %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; -[tfer_FD,~,n] = tfer_PMA.tfer_FD([],... +[tfer_FD,~,n] = tfer_pma.tfer_FD([],... m,d,1,prop,'V',V,'omega',omega); t(1) = toc; @@ -34,7 +34,7 @@ %=========================================================================% %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% -B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +B = tfer_pma.dm2zp(d,z,prop.T,prop.p); tau = B.*m; D = prop.D(B); sig = sqrt(2.*prop.L.*D./prop.v_bar); @@ -45,18 +45,18 @@ %-- Plug flow ------------------------------------------------------------% %-- Method 1S ------------------------------% tic; -[tfer_1S,G0_1S] = tfer_PMA.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); +[tfer_1S,G0_1S] = tfer_pma.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); t(2) = toc; %-- Method 1S, Ehara et al. ----------------% -tfer_Ehara = tfer_PMA.tfer_Ehara([],m,d,z,prop,'V',V,'omega',omega); +tfer_Ehara = tfer_pma.tfer_Ehara([],m,d,z,prop,'V',V,'omega',omega); %-- Parabolic flow -------------------------------------------------------% %-- Method 1S ------------------------------% tic; -[tfer_1S_pb,G0_1S_pb] = tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); +[tfer_1S_pb,G0_1S_pb] = tfer_pma.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); t(8) = toc; diff --git a/main_charge.m b/main_charge.m index 3840786..0ade0df 100644 --- a/main_charge.m +++ b/main_charge.m @@ -23,14 +23,14 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density - prop = tfer_PMA.prop_PMA('Olfert'); % get properties of the CPMA + prop = tfer_pma.prop_PMA('Olfert'); % get properties of the CPMA % prop.omega_hat = 1; % NOTE: Uncomment for APM condition %=========================================================================% %-- Finite difference solution -------------------------------------------% tic; - [tfer_FD(:,zz),sp,n{zz}] = tfer_PMA.tfer_FD(m_star,... + [tfer_FD(:,zz),sp,n{zz}] = tfer_pma.tfer_FD(m_star,... m,d,z,prop,'Rm',Rm); t(1) = toc; @@ -39,9 +39,9 @@ %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% if ~exist('d','var') - B = tfer_PMA.mp2zp(m,z,prop.T,prop.p); + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); else - B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); end tau = B.*m; D = prop.D(B); @@ -55,7 +55,7 @@ %-- Method 1C --------------------------------% tic; tfer_1C_diff(:,zz) = ... - tfer_PMA.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); + tfer_pma.tfer_1C_diff(m_star,m,d,z,prop,'Rm',Rm); t(12) = toc; end diff --git a/main_kuwata.m b/main_kuwata.m index 47b08ee..46a867d 100644 --- a/main_kuwata.m +++ b/main_kuwata.m @@ -20,14 +20,14 @@ d = 100e-9.*ones(size(m)); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_PMA('Kuwata'); % get properties of the CPMA +prop = tfer_pma.prop_PMA('Kuwata'); % get properties of the CPMA prop.D = @(B) 1e-10.*ones(size(B)); %=========================================================================% %-- Finite difference solution -------------------------------------------% -[tfer_FD,sp] = tfer_PMA.tfer_FD([],... +[tfer_FD,sp] = tfer_pma.tfer_FD([],... m,d,1,prop,'V',V,'omega',omega); @@ -35,16 +35,16 @@ %=========================================================================% %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% -B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +B = tfer_pma.dm2zp(d,z,prop.T,prop.p); tau = B.*m; %-- Particle tracking approaches -----------------------------------------% %-- Plug flow ------------------------------------------------------------% %-- Method 1S ------------------------------% [tfer_1S] = ... - tfer_PMA.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); + tfer_pma.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); [tfer_1S_pb] = ... - tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + tfer_pma.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); diff --git a/main_olfert_collings.m b/main_olfert_collings.m index 92550fd..f5685cc 100644 --- a/main_olfert_collings.m +++ b/main_olfert_collings.m @@ -21,7 +21,7 @@ d = (6.*m./(rho_eff.*pi)).^(1/3); % specify mobility diameter vector with constant effective density -prop = tfer_PMA.prop_PMA('Olfert-Collings'); % get properties of the CPMA +prop = tfer_pma.prop_PMA('Olfert-Collings'); % get properties of the CPMA prop.D = @(B) 1e-10.*ones(size(B)); omega_hat = prop.omega_hat; % only valid for CPMA @@ -30,11 +30,11 @@ %=========================================================================% %-- Finite difference solution -------------------------------------------% prop.omega_hat = 1; -[tfer_FD_w1,sp] = tfer_PMA.tfer_FD([],... +[tfer_FD_w1,sp] = tfer_pma.tfer_FD([],... m,d,1,prop,'V',V,'omega',omega); prop.omega_hat = omega_hat; -[tfer_FD,sp_cpma] = tfer_PMA.tfer_FD([],... +[tfer_FD,sp_cpma] = tfer_pma.tfer_FD([],... m,d,1,prop,'V',V,'omega',omega); @@ -42,8 +42,8 @@ %=========================================================================% %-- Transfer functions for different cases -------------------------------% %-- Setup for centriputal force ------------------------------------------% -prop = tfer_PMA.prop_CPMA('Olfert-Collings'); % get properties of the CPMA -B = tfer_PMA.dm2zp(d,z,prop.T,prop.p); +prop = tfer_pma.prop_CPMA('Olfert-Collings'); % get properties of the CPMA +B = tfer_pma.dm2zp(d,z,prop.T,prop.p); tau = B.*m; %-- Particle tracking approaches -----------------------------------------% @@ -51,11 +51,11 @@ %-- Method 1S ------------------------------% prop.omega_hat = 1; % NOTE: Uncomment for APM condition [tfer_1S_w1] = ... - tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + tfer_pma.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); prop.omega_hat = omega_hat; [tfer_1S] = ... - tfer_PMA.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); + tfer_pma.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); From 63ed7b754e4a8b638f11bb0ad6a35395f6976aeb Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 19 Sep 2019 15:22:12 -0700 Subject: [PATCH 16/22] Update to README --- README.md | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index b833a85..d7228c6 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,12 @@ -# UBC-tfer-PMA - -The attached functions and script are intended to reproduce the results of -the associated paper (submitted). They evaluate the transfer function of -the centrifugal particle mass analyzer (CPMA) and aerosol particle mass -analyzer (APM). This is done using a novel set of expressions derived from -particle tracking methods and using a finite difference method. Information -on each file is given as header information in each file, and only a brief +# mat-tfer-PMA + +The attached MATLAB functions and scripts are intended to reproduce the +results of the associated paper (submitted). They evaluate the transfer +function of particle mass analyzers (PMAs), including the centrifugal +particle mass analyzer (CPMA) and aerosol particle mass analyzer (APM). +This is done using a novel set of expressions derived from particle +tracking methods and using a finite difference method. Information on +each file is given as header information in each file, and only a brief overview is provided here. @@ -57,9 +58,10 @@ evaluates the setpoint and related properties, including C0, alpha, and beta. This script is included to demonstrate evaluation of the transfer function over multiple cases. Figure 2 that is produced by this procedure will -resemble those given in the associated work. Other scripts, `main_*.m` -are intended to replicate figures in other works and to consider multiple -charging. +resemble those given in the associated work. + +Other scripts, `main_*.m` are intended to replicate figures in other +works and to consider multiple charging. #### Remaining functions From 7399b608fdb8349291b48b19948ae754f4d381a9 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 19 Sep 2019 15:24:47 -0700 Subject: [PATCH 17/22] Forced git recognize text case --- +tfer_pma/G_fun.m | 102 ++++++++++++++++++++++++ +tfer_pma/dm2zp.m | 109 +++++++++++++++++++++++++ +tfer_pma/get_setpoint.m | 163 ++++++++++++++++++++++++++++++++++++++ +tfer_pma/mp2zp.m | 41 ++++++++++ +tfer_pma/prop_PMA.m | 92 ++++++++++++++++++++++ +tfer_pma/tfer_1C.m | 41 ++++++++++ +tfer_pma/tfer_1C_diff.m | 55 +++++++++++++ +tfer_pma/tfer_1C_pb.m | 59 ++++++++++++++ +tfer_pma/tfer_1S.m | 49 ++++++++++++ +tfer_pma/tfer_1S_diff.m | 55 +++++++++++++ +tfer_pma/tfer_1S_pb.m | 58 ++++++++++++++ +tfer_pma/tfer_2C.m | 47 +++++++++++ +tfer_pma/tfer_2C_diff.m | 55 +++++++++++++ +tfer_pma/tfer_2S.m | 54 +++++++++++++ +tfer_pma/tfer_2S_diff.m | 55 +++++++++++++ +tfer_pma/tfer_Ehara.m | 51 ++++++++++++ +tfer_pma/tfer_FD.m | 166 +++++++++++++++++++++++++++++++++++++++ +tfer_pma/tfer_GE.m | 58 ++++++++++++++ +tfer_pma/tfer_GE_diff.m | 55 +++++++++++++ +tfer_pma/tfer_W1.m | 49 ++++++++++++ +tfer_pma/tfer_W1_diff.m | 55 +++++++++++++ +tfer_pma/tfer_W1_pb.m | 59 ++++++++++++++ +tfer_pma/tfer_tri.m | 51 ++++++++++++ +tfer_pma/tridiag.m | 56 +++++++++++++ main_ehara.m | 81 +++++++++++++++++++ 25 files changed, 1716 insertions(+) create mode 100644 +tfer_pma/G_fun.m create mode 100644 +tfer_pma/dm2zp.m create mode 100644 +tfer_pma/get_setpoint.m create mode 100644 +tfer_pma/mp2zp.m create mode 100644 +tfer_pma/prop_PMA.m create mode 100644 +tfer_pma/tfer_1C.m create mode 100644 +tfer_pma/tfer_1C_diff.m create mode 100644 +tfer_pma/tfer_1C_pb.m create mode 100644 +tfer_pma/tfer_1S.m create mode 100644 +tfer_pma/tfer_1S_diff.m create mode 100644 +tfer_pma/tfer_1S_pb.m create mode 100644 +tfer_pma/tfer_2C.m create mode 100644 +tfer_pma/tfer_2C_diff.m create mode 100644 +tfer_pma/tfer_2S.m create mode 100644 +tfer_pma/tfer_2S_diff.m create mode 100644 +tfer_pma/tfer_Ehara.m create mode 100644 +tfer_pma/tfer_FD.m create mode 100644 +tfer_pma/tfer_GE.m create mode 100644 +tfer_pma/tfer_GE_diff.m create mode 100644 +tfer_pma/tfer_W1.m create mode 100644 +tfer_pma/tfer_W1_diff.m create mode 100644 +tfer_pma/tfer_W1_pb.m create mode 100644 +tfer_pma/tfer_tri.m create mode 100644 +tfer_pma/tridiag.m create mode 100644 main_ehara.m diff --git a/+tfer_pma/G_fun.m b/+tfer_pma/G_fun.m new file mode 100644 index 0000000..7bc7d94 --- /dev/null +++ b/+tfer_pma/G_fun.m @@ -0,0 +1,102 @@ + +% G_FUN Function used in root-finding procedure to specify search interval. +% Author: Timothy Sipkens, 2019-02-02 +%=========================================================================% + +function G = G_fun(min_fun,rL,rs,r1,r2,alpha,beta) +%-------------------------------------------------------------------------% +% Inputs: +% min_fun Function to minimize, depends on case letter (Sipkens et al.) +% rL Particle exit radius +% rs Equilibirum radius +% r1 Radius of inner electrode +% r2 Radius of outer electrode +% alpha First parameter from azimuthal velocity distribution +% beta Second parameter from azimuthal velocity distribution +% +% Outputs: +% G The value of G0 following root-finding procedure +%-------------------------------------------------------------------------% + + +condit = (alpha^2) < (beta^2./(rs.^4)); + % whether or not lambda is positive of negative + +ns = length(rs); % number of equilirbrium radii to consider (one per particle mass considered) +nL = length(rL); % number of final radial positions to consider (usually only one) +offset = 1e-9; % offset of rs in which rL is considered equal to rs + % (.e. no change in radial position of particle) + +G = zeros(nL,ns); % pre-allocate output variable + + +%-- Determine interval in which to search and fo root-finding ------------% +for jj=1:nL % loop throught all specified particle exit radii + for ii=1:ns % loop to consider multiple equilbirium radii (multiple m) + + if rL(jj)r2 + % particles cannot exist here, return r2 + G(jj,ii) = r2; + + elseif and(rL(jj)<(rs(ii)+5*offset),rL(jj)>(rs(ii)-5*offset)) + % if rL = rs, within offset, then r0 = rs + G(jj,ii) = rL(jj); + + elseif rL(jj)>rs(ii) + % if rL > rs, then r0 > rL + + if condit(ii) % zero is in interval [rL,r2] + if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),r2,ii)) + % if sign does not change in interval [rL,r2], return r2 + G(jj,ii) = r2; + + else % find zero in interval [rL,r2] + G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [rL(jj),r2]); + end + + else % zero is in interval [max(r1,rs),rL] + if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),max(r1,rs(ii)+offset),ii)) + % if sign does not change in interval [max(r1,rs),rL], return max(r1,rs) + G(jj,ii) = max(r1,rs(ii)+offset); + + else % find zero in interval [max(r1,rs),rL] + G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [max(r1,rs(ii)+offset),rL(jj)]); + + end + end + + + elseif rL(jj)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_pma/tfer_1C_pb.m b/+tfer_pma/tfer_1C_pb.m new file mode 100644 index 0000000..233c30b --- /dev/null +++ b/+tfer_pma/tfer_1C_pb.m @@ -0,0 +1,59 @@ + +% TFER_1C_PB Evaluates the transfer function for a PMA in Case B (w/ parabolic flow). +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_1C_pb(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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); + + +%-- 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 + + +%-- 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) tfer_pma.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; + +end + diff --git a/+tfer_pma/tfer_1S.m b/+tfer_pma/tfer_1S.m new file mode 100644 index 0000000..caaa97e --- /dev/null +++ b/+tfer_pma/tfer_1S.m @@ -0,0 +1,49 @@ + +% TFER_1S Evaluates the transfer function for a PMA in Case A. +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_1S(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 % else, use other root + rs = real((sqrt(C0./m)+sqrt(C0./m-4*sp.alpha*sp.beta))./(2*sp.alpha)); +end + + +%-- Estimate device parameter --------------------------------------------% +lam = 2.*tau.*(sp.alpha^2-sp.beta^2./(rs.^4)).*prop.L./prop.v_bar; + + +%-- Evaluate G0 and transfer function ------------------------------------% +G0 = @(r) rs+(r-rs).*exp(-lam); + +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); + +end + diff --git a/+tfer_pma/tfer_1S_diff.m b/+tfer_pma/tfer_1S_diff.m new file mode 100644 index 0000000..fadcb96 --- /dev/null +++ b/+tfer_pma/tfer_1S_diff.m @@ -0,0 +1,55 @@ + +% TFER_1S_DIFF Evaluates the transfer function for a PMA in Case A (w/ diffusion). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_1S_diff(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +%-- Evaluate mechanical mobility for diffusion calc. ---------------------% +if ~exist('d','var') + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); + % if mobility is not specified, use mass-mobility relation to estimate +else + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); +end + +D = prop.D(B).*z; + % diffusion coefficient is previously defined function multiplied by + % integer charge state +sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter + +[~,G0] = tfer_pma.tfer_1S(m_star,m,d,z,prop,varargin{:}); + % get G0 function for this case + +rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 + +end diff --git a/+tfer_pma/tfer_1S_pb.m b/+tfer_pma/tfer_1S_pb.m new file mode 100644 index 0000000..d0eefdf --- /dev/null +++ b/+tfer_pma/tfer_1S_pb.m @@ -0,0 +1,58 @@ + +% TFER_1S_PB Evaluates the transfer function for a PMA in Case A (w/ parabolic flow). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_1S_pb(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 device parameter --------------------------------------------% +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; + + +%-- 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) tfer_pma.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; + +end + diff --git a/+tfer_pma/tfer_2C.m b/+tfer_pma/tfer_2C.m new file mode 100644 index 0000000..0b75275 --- /dev/null +++ b/+tfer_pma/tfer_2C.m @@ -0,0 +1,47 @@ + +% TFER_2C Evaluates the transfer function for a PMA in Case D. +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_2C(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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))); + +C6 = -1./(sqrt(4.*C3.*C5-C4.^2)); + +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; + +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)); + +end + diff --git a/+tfer_pma/tfer_2C_diff.m b/+tfer_pma/tfer_2C_diff.m new file mode 100644 index 0000000..97120e3 --- /dev/null +++ b/+tfer_pma/tfer_2C_diff.m @@ -0,0 +1,55 @@ + +% TFER_2C_DIFF Evaluates the transfer function for a PMA in Case D (w/ diffusion). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_2C_diff(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +%-- Evaluate mechanical mobility for diffusion calc. ---------------------% +if ~exist('d','var') + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); + % if mobility is not specified, use mass-mobility relation to estimate +else + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); +end + +D = prop.D(B).*z; + % diffusion coefficient is previously defined function multiplied by + % integer charge state +sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter + +[~,G0] = tfer_pma.tfer_2C(m_star,m,d,z,prop,varargin{:}); + % get G0 function for this case + +rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 + +end diff --git a/+tfer_pma/tfer_2S.m b/+tfer_pma/tfer_2S.m new file mode 100644 index 0000000..06ecc44 --- /dev/null +++ b/+tfer_pma/tfer_2S.m @@ -0,0 +1,54 @@ + +% TFER_2S Evaluates the transfer function for a PMA in Case C. +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_2S(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 device parameter --------------------------------------------% +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)); + + +%-- Evaluate G0 and transfer function ------------------------------------% +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))); + +Lambda = (1/(2*prop.del)).*(rb-ra); + +end + diff --git a/+tfer_pma/tfer_2S_diff.m b/+tfer_pma/tfer_2S_diff.m new file mode 100644 index 0000000..9045298 --- /dev/null +++ b/+tfer_pma/tfer_2S_diff.m @@ -0,0 +1,55 @@ + +% TFER_2S_DIFF Evaluates the transfer function for a PMA in Case C (w/ diffusion). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_2S_diff(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +%-- Evaluate mechanical mobility for diffusion calc. ---------------------% +if ~exist('d','var') + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); + % if mobility is not specified, use mass-mobility relation to estimate +else + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); +end + +D = prop.D(B).*z; + % diffusion coefficient is previously defined function multiplied by + % integer charge state +sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter + +[~,G0] = tfer_pma.tfer_2S(m_star,m,d,z,prop,varargin{:}); + % get G0 function for this case + +rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 + +end diff --git a/+tfer_pma/tfer_Ehara.m b/+tfer_pma/tfer_Ehara.m new file mode 100644 index 0000000..f880d71 --- /dev/null +++ b/+tfer_pma/tfer_Ehara.m @@ -0,0 +1,51 @@ + +% TFER_EHARA Evaluates the transfer function for a PMA in Case A as per Ehara et al. (1996). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda] = tfer_Ehara(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 device parameter --------------------------------------------% +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(1(z.*m_star+2.*sp.m_max),... + m<(z.*m_star-2.*sp.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)); + + + %-- 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 nargout==3 % 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 = tfer_pma.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==3; n_mat(jj,:) = n_vec; end + % record particle distribution at this axial position + end + + %-- Format output ------------------------------------------% + if nargout==3 % for visualizing number concentrations + n.n_mat{ii} = max(n_mat,0); + 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 + + +%== 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) + +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)]; + + +%-- 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 + +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_pma/tfer_GE.m b/+tfer_pma/tfer_GE.m new file mode 100644 index 0000000..9d0b3fe --- /dev/null +++ b/+tfer_pma/tfer_GE.m @@ -0,0 +1,58 @@ + +% TFER_GE Evaluates the transfer function for a PMA in Case F. +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_GE(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 ----------------------------------------% +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; + + +%-- 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; + + +%-- Evaluate G0 and transfer function ------------------------------------% +G0 = @(r) tfer_pma.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 = (1/(2*prop.del)).*(rb-ra); + +end + diff --git a/+tfer_pma/tfer_GE_diff.m b/+tfer_pma/tfer_GE_diff.m new file mode 100644 index 0000000..015974b --- /dev/null +++ b/+tfer_pma/tfer_GE_diff.m @@ -0,0 +1,55 @@ + +% TFER_GE_DIFF Evaluates the transfer function for a PMA in Case F (w/ diffusion). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_GE_diff(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +%-- Evaluate mechanical mobility for diffusion calc. ---------------------% +if ~exist('d','var') + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); + % if mobility is not specified, use mass-mobility relation to estimate +else + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); +end + +D = prop.D(B).*z; + % diffusion coefficient is previously defined function multiplied by + % integer charge state +sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter + +[~,G0] = tfer_pma.tfer_GE(m_star,m,d,z,prop,varargin{:}); + % get G0 function for this case + +rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 + +end diff --git a/+tfer_pma/tfer_W1.m b/+tfer_pma/tfer_W1.m new file mode 100644 index 0000000..b9751a7 --- /dev/null +++ b/+tfer_pma/tfer_W1.m @@ -0,0 +1,49 @@ + +% TFER_W1 Evaluates the transfer function for a PMA in Case E. +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_W1(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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 device parameter --------------------------------------------% +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); + +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); + +end + diff --git a/+tfer_pma/tfer_W1_diff.m b/+tfer_pma/tfer_W1_diff.m new file mode 100644 index 0000000..d3741eb --- /dev/null +++ b/+tfer_pma/tfer_W1_diff.m @@ -0,0 +1,55 @@ + +% TFER_W1_DIFF Evaluates the transfer function for a PMA in Case E (w/ diffusion). +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda,G0] = tfer_W1_diff(m_star,m,d,z,prop,varargin) +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +%-- Evaluate mechanical mobility for diffusion calc. ---------------------% +if ~exist('d','var') + B = tfer_pma.mp2zp(m,z,prop.T,prop.p); + % if mobility is not specified, use mass-mobility relation to estimate +else + B = tfer_pma.dm2zp(d,z,prop.T,prop.p); +end + +D = prop.D(B).*z; + % diffusion coefficient is previously defined function multiplied by + % integer charge state +sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter + +[~,G0] = tfer_pma.tfer_W1(m_star,m,d,z,prop,varargin{:}); + % get G0 function for this case + +rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 + +end diff --git a/+tfer_pma/tfer_W1_pb.m b/+tfer_pma/tfer_W1_pb.m new file mode 100644 index 0000000..b84ecc6 --- /dev/null +++ b/+tfer_pma/tfer_W1_pb.m @@ -0,0 +1,59 @@ + +% TFER_W1_PB Evaluates the transfer function for a PMA in Case E (w/ parabolic flow). +% Author: Timothy Sipkens, 2019-03-21 +%=========================================================================% + +function [Lambda,G0] = tfer_W1_pb(m_star,m,d,z,prop,varargin) +% +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Device properties (e.g. classifier length) +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Outputs: +% Lambda Transfer function +% G0 Function mapping final to initial radial position +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint (parses d and z) + +%-- 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) tfer_pma.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; + +end + diff --git a/+tfer_pma/tfer_tri.m b/+tfer_pma/tfer_tri.m new file mode 100644 index 0000000..69ca2da --- /dev/null +++ b/+tfer_pma/tfer_tri.m @@ -0,0 +1,51 @@ + +% TFER_TRI Evaluates the transfer function for a PMA as a triangular function. +% Author: Timothy Sipkens, 2018-12-27 +%=========================================================================% + +function [Lambda] = tfer_tri(m_star,m,d,z,prop,varargin) +% +%-------------------------------------------------------------------------% +% Inputs: +% m_star Setpoint particle mass +% m Particle mass +% d Particle mobility diameter +% z Integer charge state +% prop Properties of particle mass analyzer +% varargin Name-value pairs for setpoint (Optional, default Rm = 3) +% ('Rm',double) - Resolution +% ('omega1',double) - Angular speed of inner electrode +% ('V',double) - Setpoint voltage +% +% Output: +% Lambda Transfer function +%-------------------------------------------------------------------------% + + +tfer_pma.get_setpoint; % get setpoint + +if ~isfield(sp,'m_max') % if m_max was not specified + n_B = -0.6436; + B_star = mp2zp(m_star,1,prop.T,prop.p); % involves invoking mass-mobility relation + + omega = sp.omega1.*... + ((prop.r_hat^2-prop.omega_hat)/(prop.r_hat^2-1)+... + prop.r1^2*(prop.omega_hat-1)/(prop.r_hat^2-1)/prop.rc^2); + + m_ratio = fzero(@(x) abs(x).^(n_B+1)-abs(x).^n_B-... + prop.Q/(omega.^2.*m_star*B_star*2*pi*prop.rc^2*prop.L),m_star.*1.5); + sp.m_max = m_star.*abs(m_ratio); + +end + + +m_del = sp.m_max-m_star; % FWHM of the transfer function (related to resolution) +m_min = 2.*m_star-sp.m_max; % lower end of the transfer function + +Lambda = zeros(size(m))+... + (m<=m_star).*(m>m_min).*(m-m_min)./m_del+... + (m>m_star).*(m Date: Thu, 19 Sep 2019 15:56:34 -0700 Subject: [PATCH 18/22] Resolving duplicate folder names (step 1) --- +tfer_pma/G_fun.m | 102 -------------- +tfer_pma/dm2zp.m | 109 --------------- +tfer_pma/get_setpoint.m | 163 ---------------------- +tfer_pma/mp2zp.m | 41 ------ +tfer_pma/prop_PMA.m | 92 ------------- +tfer_pma/tfer_1C.m | 41 ------ +tfer_pma/tfer_1C_diff.m | 55 -------- +tfer_pma/tfer_1C_pb.m | 59 -------- +tfer_pma/tfer_1S.m | 49 ------- +tfer_pma/tfer_1S_diff.m | 55 -------- +tfer_pma/tfer_1S_pb.m | 58 -------- +tfer_pma/tfer_2C.m | 47 ------- +tfer_pma/tfer_2C_diff.m | 55 -------- +tfer_pma/tfer_2S.m | 54 -------- +tfer_pma/tfer_2S_diff.m | 55 -------- +tfer_pma/tfer_Ehara.m | 51 ------- +tfer_pma/tfer_FD.m | 166 ----------------------- +tfer_pma/tfer_GE.m | 58 -------- +tfer_pma/tfer_GE_diff.m | 55 -------- +tfer_pma/tfer_W1.m | 49 ------- +tfer_pma/tfer_W1_diff.m | 55 -------- +tfer_pma/tfer_W1_pb.m | 59 -------- +tfer_pma/tfer_tri.m | 51 ------- +tfer_pma/tridiag.m | 56 -------- {+tfer_PMA => +tfer_pma2}/G_fun.m | 0 {+tfer_PMA => +tfer_pma2}/dm2zp.m | 0 {+tfer_PMA => +tfer_pma2}/get_setpoint.m | 0 {+tfer_PMA => +tfer_pma2}/mp2zp.m | 0 {+tfer_PMA => +tfer_pma2}/prop_PMA.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1C.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1C_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1C_pb.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1S.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1S_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_1S_pb.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_2C.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_2C_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_2S.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_2S_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_Ehara.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_FD.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_GE.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_GE_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_W1.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_W1_diff.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_W1_pb.m | 0 {+tfer_PMA => +tfer_pma2}/tfer_tri.m | 0 {+tfer_PMA => +tfer_pma2}/tridiag.m | 0 48 files changed, 1635 deletions(-) delete mode 100644 +tfer_pma/G_fun.m delete mode 100644 +tfer_pma/dm2zp.m delete mode 100644 +tfer_pma/get_setpoint.m delete mode 100644 +tfer_pma/mp2zp.m delete mode 100644 +tfer_pma/prop_PMA.m delete mode 100644 +tfer_pma/tfer_1C.m delete mode 100644 +tfer_pma/tfer_1C_diff.m delete mode 100644 +tfer_pma/tfer_1C_pb.m delete mode 100644 +tfer_pma/tfer_1S.m delete mode 100644 +tfer_pma/tfer_1S_diff.m delete mode 100644 +tfer_pma/tfer_1S_pb.m delete mode 100644 +tfer_pma/tfer_2C.m delete mode 100644 +tfer_pma/tfer_2C_diff.m delete mode 100644 +tfer_pma/tfer_2S.m delete mode 100644 +tfer_pma/tfer_2S_diff.m delete mode 100644 +tfer_pma/tfer_Ehara.m delete mode 100644 +tfer_pma/tfer_FD.m delete mode 100644 +tfer_pma/tfer_GE.m delete mode 100644 +tfer_pma/tfer_GE_diff.m delete mode 100644 +tfer_pma/tfer_W1.m delete mode 100644 +tfer_pma/tfer_W1_diff.m delete mode 100644 +tfer_pma/tfer_W1_pb.m delete mode 100644 +tfer_pma/tfer_tri.m delete mode 100644 +tfer_pma/tridiag.m rename {+tfer_PMA => +tfer_pma2}/G_fun.m (100%) rename {+tfer_PMA => +tfer_pma2}/dm2zp.m (100%) rename {+tfer_PMA => +tfer_pma2}/get_setpoint.m (100%) rename {+tfer_PMA => +tfer_pma2}/mp2zp.m (100%) rename {+tfer_PMA => +tfer_pma2}/prop_PMA.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1C.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1C_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1C_pb.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1S.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1S_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_1S_pb.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_2C.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_2C_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_2S.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_2S_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_Ehara.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_FD.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_GE.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_GE_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_W1.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_W1_diff.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_W1_pb.m (100%) rename {+tfer_PMA => +tfer_pma2}/tfer_tri.m (100%) rename {+tfer_PMA => +tfer_pma2}/tridiag.m (100%) diff --git a/+tfer_pma/G_fun.m b/+tfer_pma/G_fun.m deleted file mode 100644 index 7bc7d94..0000000 --- a/+tfer_pma/G_fun.m +++ /dev/null @@ -1,102 +0,0 @@ - -% G_FUN Function used in root-finding procedure to specify search interval. -% Author: Timothy Sipkens, 2019-02-02 -%=========================================================================% - -function G = G_fun(min_fun,rL,rs,r1,r2,alpha,beta) -%-------------------------------------------------------------------------% -% Inputs: -% min_fun Function to minimize, depends on case letter (Sipkens et al.) -% rL Particle exit radius -% rs Equilibirum radius -% r1 Radius of inner electrode -% r2 Radius of outer electrode -% alpha First parameter from azimuthal velocity distribution -% beta Second parameter from azimuthal velocity distribution -% -% Outputs: -% G The value of G0 following root-finding procedure -%-------------------------------------------------------------------------% - - -condit = (alpha^2) < (beta^2./(rs.^4)); - % whether or not lambda is positive of negative - -ns = length(rs); % number of equilirbrium radii to consider (one per particle mass considered) -nL = length(rL); % number of final radial positions to consider (usually only one) -offset = 1e-9; % offset of rs in which rL is considered equal to rs - % (.e. no change in radial position of particle) - -G = zeros(nL,ns); % pre-allocate output variable - - -%-- Determine interval in which to search and fo root-finding ------------% -for jj=1:nL % loop throught all specified particle exit radii - for ii=1:ns % loop to consider multiple equilbirium radii (multiple m) - - if rL(jj)r2 - % particles cannot exist here, return r2 - G(jj,ii) = r2; - - elseif and(rL(jj)<(rs(ii)+5*offset),rL(jj)>(rs(ii)-5*offset)) - % if rL = rs, within offset, then r0 = rs - G(jj,ii) = rL(jj); - - elseif rL(jj)>rs(ii) - % if rL > rs, then r0 > rL - - if condit(ii) % zero is in interval [rL,r2] - if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),r2,ii)) - % if sign does not change in interval [rL,r2], return r2 - G(jj,ii) = r2; - - else % find zero in interval [rL,r2] - G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [rL(jj),r2]); - end - - else % zero is in interval [max(r1,rs),rL] - if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),max(r1,rs(ii)+offset),ii)) - % if sign does not change in interval [max(r1,rs),rL], return max(r1,rs) - G(jj,ii) = max(r1,rs(ii)+offset); - - else % find zero in interval [max(r1,rs),rL] - G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [max(r1,rs(ii)+offset),rL(jj)]); - - end - end - - - elseif rL(jj)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_pma/tfer_1C_pb.m b/+tfer_pma/tfer_1C_pb.m deleted file mode 100644 index 233c30b..0000000 --- a/+tfer_pma/tfer_1C_pb.m +++ /dev/null @@ -1,59 +0,0 @@ - -% TFER_1C_PB Evaluates the transfer function for a PMA in Case B (w/ parabolic flow). -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_1C_pb(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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); - - -%-- 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 - - -%-- 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) tfer_pma.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; - -end - diff --git a/+tfer_pma/tfer_1S.m b/+tfer_pma/tfer_1S.m deleted file mode 100644 index caaa97e..0000000 --- a/+tfer_pma/tfer_1S.m +++ /dev/null @@ -1,49 +0,0 @@ - -% TFER_1S Evaluates the transfer function for a PMA in Case A. -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_1S(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 % else, use other root - rs = real((sqrt(C0./m)+sqrt(C0./m-4*sp.alpha*sp.beta))./(2*sp.alpha)); -end - - -%-- Estimate device parameter --------------------------------------------% -lam = 2.*tau.*(sp.alpha^2-sp.beta^2./(rs.^4)).*prop.L./prop.v_bar; - - -%-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) rs+(r-rs).*exp(-lam); - -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); - -end - diff --git a/+tfer_pma/tfer_1S_diff.m b/+tfer_pma/tfer_1S_diff.m deleted file mode 100644 index fadcb96..0000000 --- a/+tfer_pma/tfer_1S_diff.m +++ /dev/null @@ -1,55 +0,0 @@ - -% TFER_1S_DIFF Evaluates the transfer function for a PMA in Case A (w/ diffusion). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_1S_diff(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -%-- Evaluate mechanical mobility for diffusion calc. ---------------------% -if ~exist('d','var') - B = tfer_pma.mp2zp(m,z,prop.T,prop.p); - % if mobility is not specified, use mass-mobility relation to estimate -else - B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -end - -D = prop.D(B).*z; - % diffusion coefficient is previously defined function multiplied by - % integer charge state -sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter - -[~,G0] = tfer_pma.tfer_1S(m_star,m,d,z,prop,varargin{:}); - % get G0 function for this case - -rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 - -end diff --git a/+tfer_pma/tfer_1S_pb.m b/+tfer_pma/tfer_1S_pb.m deleted file mode 100644 index d0eefdf..0000000 --- a/+tfer_pma/tfer_1S_pb.m +++ /dev/null @@ -1,58 +0,0 @@ - -% TFER_1S_PB Evaluates the transfer function for a PMA in Case A (w/ parabolic flow). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_1S_pb(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 device parameter --------------------------------------------% -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; - - -%-- 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) tfer_pma.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; - -end - diff --git a/+tfer_pma/tfer_2C.m b/+tfer_pma/tfer_2C.m deleted file mode 100644 index 0b75275..0000000 --- a/+tfer_pma/tfer_2C.m +++ /dev/null @@ -1,47 +0,0 @@ - -% TFER_2C Evaluates the transfer function for a PMA in Case D. -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_2C(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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))); - -C6 = -1./(sqrt(4.*C3.*C5-C4.^2)); - -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; - -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)); - -end - diff --git a/+tfer_pma/tfer_2C_diff.m b/+tfer_pma/tfer_2C_diff.m deleted file mode 100644 index 97120e3..0000000 --- a/+tfer_pma/tfer_2C_diff.m +++ /dev/null @@ -1,55 +0,0 @@ - -% TFER_2C_DIFF Evaluates the transfer function for a PMA in Case D (w/ diffusion). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_2C_diff(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -%-- Evaluate mechanical mobility for diffusion calc. ---------------------% -if ~exist('d','var') - B = tfer_pma.mp2zp(m,z,prop.T,prop.p); - % if mobility is not specified, use mass-mobility relation to estimate -else - B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -end - -D = prop.D(B).*z; - % diffusion coefficient is previously defined function multiplied by - % integer charge state -sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter - -[~,G0] = tfer_pma.tfer_2C(m_star,m,d,z,prop,varargin{:}); - % get G0 function for this case - -rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 - -end diff --git a/+tfer_pma/tfer_2S.m b/+tfer_pma/tfer_2S.m deleted file mode 100644 index 06ecc44..0000000 --- a/+tfer_pma/tfer_2S.m +++ /dev/null @@ -1,54 +0,0 @@ - -% TFER_2S Evaluates the transfer function for a PMA in Case C. -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_2S(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 device parameter --------------------------------------------% -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)); - - -%-- Evaluate G0 and transfer function ------------------------------------% -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))); - -Lambda = (1/(2*prop.del)).*(rb-ra); - -end - diff --git a/+tfer_pma/tfer_2S_diff.m b/+tfer_pma/tfer_2S_diff.m deleted file mode 100644 index 9045298..0000000 --- a/+tfer_pma/tfer_2S_diff.m +++ /dev/null @@ -1,55 +0,0 @@ - -% TFER_2S_DIFF Evaluates the transfer function for a PMA in Case C (w/ diffusion). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_2S_diff(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -%-- Evaluate mechanical mobility for diffusion calc. ---------------------% -if ~exist('d','var') - B = tfer_pma.mp2zp(m,z,prop.T,prop.p); - % if mobility is not specified, use mass-mobility relation to estimate -else - B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -end - -D = prop.D(B).*z; - % diffusion coefficient is previously defined function multiplied by - % integer charge state -sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter - -[~,G0] = tfer_pma.tfer_2S(m_star,m,d,z,prop,varargin{:}); - % get G0 function for this case - -rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 - -end diff --git a/+tfer_pma/tfer_Ehara.m b/+tfer_pma/tfer_Ehara.m deleted file mode 100644 index f880d71..0000000 --- a/+tfer_pma/tfer_Ehara.m +++ /dev/null @@ -1,51 +0,0 @@ - -% TFER_EHARA Evaluates the transfer function for a PMA in Case A as per Ehara et al. (1996). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda] = tfer_Ehara(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 device parameter --------------------------------------------% -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(1(z.*m_star+2.*sp.m_max),... - m<(z.*m_star-2.*sp.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)); - - - %-- 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 nargout==3 % 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 = tfer_pma.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==3; n_mat(jj,:) = n_vec; end - % record particle distribution at this axial position - end - - %-- Format output ------------------------------------------% - if nargout==3 % for visualizing number concentrations - n.n_mat{ii} = max(n_mat,0); - 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 - - -%== 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) - -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)]; - - -%-- 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 - -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_pma/tfer_GE.m b/+tfer_pma/tfer_GE.m deleted file mode 100644 index 9d0b3fe..0000000 --- a/+tfer_pma/tfer_GE.m +++ /dev/null @@ -1,58 +0,0 @@ - -% TFER_GE Evaluates the transfer function for a PMA in Case F. -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_GE(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 ----------------------------------------% -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; - - -%-- 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; - - -%-- Evaluate G0 and transfer function ------------------------------------% -G0 = @(r) tfer_pma.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 = (1/(2*prop.del)).*(rb-ra); - -end - diff --git a/+tfer_pma/tfer_GE_diff.m b/+tfer_pma/tfer_GE_diff.m deleted file mode 100644 index 015974b..0000000 --- a/+tfer_pma/tfer_GE_diff.m +++ /dev/null @@ -1,55 +0,0 @@ - -% TFER_GE_DIFF Evaluates the transfer function for a PMA in Case F (w/ diffusion). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_GE_diff(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -%-- Evaluate mechanical mobility for diffusion calc. ---------------------% -if ~exist('d','var') - B = tfer_pma.mp2zp(m,z,prop.T,prop.p); - % if mobility is not specified, use mass-mobility relation to estimate -else - B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -end - -D = prop.D(B).*z; - % diffusion coefficient is previously defined function multiplied by - % integer charge state -sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter - -[~,G0] = tfer_pma.tfer_GE(m_star,m,d,z,prop,varargin{:}); - % get G0 function for this case - -rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 - -end diff --git a/+tfer_pma/tfer_W1.m b/+tfer_pma/tfer_W1.m deleted file mode 100644 index b9751a7..0000000 --- a/+tfer_pma/tfer_W1.m +++ /dev/null @@ -1,49 +0,0 @@ - -% TFER_W1 Evaluates the transfer function for a PMA in Case E. -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_W1(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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 device parameter --------------------------------------------% -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); - -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); - -end - diff --git a/+tfer_pma/tfer_W1_diff.m b/+tfer_pma/tfer_W1_diff.m deleted file mode 100644 index d3741eb..0000000 --- a/+tfer_pma/tfer_W1_diff.m +++ /dev/null @@ -1,55 +0,0 @@ - -% TFER_W1_DIFF Evaluates the transfer function for a PMA in Case E (w/ diffusion). -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda,G0] = tfer_W1_diff(m_star,m,d,z,prop,varargin) -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -%-- Evaluate mechanical mobility for diffusion calc. ---------------------% -if ~exist('d','var') - B = tfer_pma.mp2zp(m,z,prop.T,prop.p); - % if mobility is not specified, use mass-mobility relation to estimate -else - B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -end - -D = prop.D(B).*z; - % diffusion coefficient is previously defined function multiplied by - % integer charge state -sig = sqrt(2.*prop.L.*D./prop.v_bar); % diffusive spreading parameter - -[~,G0] = tfer_pma.tfer_W1(m_star,m,d,z,prop,varargin{:}); - % get G0 function for this case - -rho_fun = @(G,r) (G-r)./(sqrt(2).*sig); % reuccring 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 - -end diff --git a/+tfer_pma/tfer_W1_pb.m b/+tfer_pma/tfer_W1_pb.m deleted file mode 100644 index b84ecc6..0000000 --- a/+tfer_pma/tfer_W1_pb.m +++ /dev/null @@ -1,59 +0,0 @@ - -% TFER_W1_PB Evaluates the transfer function for a PMA in Case E (w/ parabolic flow). -% Author: Timothy Sipkens, 2019-03-21 -%=========================================================================% - -function [Lambda,G0] = tfer_W1_pb(m_star,m,d,z,prop,varargin) -% -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Device properties (e.g. classifier length) -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Outputs: -% Lambda Transfer function -% G0 Function mapping final to initial radial position -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint (parses d and z) - -%-- 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) tfer_pma.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; - -end - diff --git a/+tfer_pma/tfer_tri.m b/+tfer_pma/tfer_tri.m deleted file mode 100644 index 69ca2da..0000000 --- a/+tfer_pma/tfer_tri.m +++ /dev/null @@ -1,51 +0,0 @@ - -% TFER_TRI Evaluates the transfer function for a PMA as a triangular function. -% Author: Timothy Sipkens, 2018-12-27 -%=========================================================================% - -function [Lambda] = tfer_tri(m_star,m,d,z,prop,varargin) -% -%-------------------------------------------------------------------------% -% Inputs: -% m_star Setpoint particle mass -% m Particle mass -% d Particle mobility diameter -% z Integer charge state -% prop Properties of particle mass analyzer -% varargin Name-value pairs for setpoint (Optional, default Rm = 3) -% ('Rm',double) - Resolution -% ('omega1',double) - Angular speed of inner electrode -% ('V',double) - Setpoint voltage -% -% Output: -% Lambda Transfer function -%-------------------------------------------------------------------------% - - -tfer_pma.get_setpoint; % get setpoint - -if ~isfield(sp,'m_max') % if m_max was not specified - n_B = -0.6436; - B_star = mp2zp(m_star,1,prop.T,prop.p); % involves invoking mass-mobility relation - - omega = sp.omega1.*... - ((prop.r_hat^2-prop.omega_hat)/(prop.r_hat^2-1)+... - prop.r1^2*(prop.omega_hat-1)/(prop.r_hat^2-1)/prop.rc^2); - - m_ratio = fzero(@(x) abs(x).^(n_B+1)-abs(x).^n_B-... - prop.Q/(omega.^2.*m_star*B_star*2*pi*prop.rc^2*prop.L),m_star.*1.5); - sp.m_max = m_star.*abs(m_ratio); - -end - - -m_del = sp.m_max-m_star; % FWHM of the transfer function (related to resolution) -m_min = 2.*m_star-sp.m_max; % lower end of the transfer function - -Lambda = zeros(size(m))+... - (m<=m_star).*(m>m_min).*(m-m_min)./m_del+... - (m>m_star).*(m Date: Thu, 19 Sep 2019 15:57:41 -0700 Subject: [PATCH 19/22] Revert back to original naming --- {+tfer_pma2 => +tfer_pma}/G_fun.m | 0 {+tfer_pma2 => +tfer_pma}/dm2zp.m | 0 {+tfer_pma2 => +tfer_pma}/get_setpoint.m | 0 {+tfer_pma2 => +tfer_pma}/mp2zp.m | 0 {+tfer_pma2 => +tfer_pma}/prop_PMA.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1C.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1C_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1C_pb.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1S.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1S_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_1S_pb.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_2C.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_2C_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_2S.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_2S_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_Ehara.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_FD.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_GE.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_GE_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_W1.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_W1_diff.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_W1_pb.m | 0 {+tfer_pma2 => +tfer_pma}/tfer_tri.m | 0 {+tfer_pma2 => +tfer_pma}/tridiag.m | 0 main_ehara.m | 81 ------------------------ main_Ehara.m => main_ehara2.m | 0 26 files changed, 81 deletions(-) rename {+tfer_pma2 => +tfer_pma}/G_fun.m (100%) rename {+tfer_pma2 => +tfer_pma}/dm2zp.m (100%) rename {+tfer_pma2 => +tfer_pma}/get_setpoint.m (100%) rename {+tfer_pma2 => +tfer_pma}/mp2zp.m (100%) rename {+tfer_pma2 => +tfer_pma}/prop_PMA.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1C.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1C_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1C_pb.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1S.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1S_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_1S_pb.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_2C.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_2C_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_2S.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_2S_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_Ehara.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_FD.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_GE.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_GE_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_W1.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_W1_diff.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_W1_pb.m (100%) rename {+tfer_pma2 => +tfer_pma}/tfer_tri.m (100%) rename {+tfer_pma2 => +tfer_pma}/tridiag.m (100%) delete mode 100644 main_ehara.m rename main_Ehara.m => main_ehara2.m (100%) diff --git a/+tfer_pma2/G_fun.m b/+tfer_pma/G_fun.m similarity index 100% rename from +tfer_pma2/G_fun.m rename to +tfer_pma/G_fun.m diff --git a/+tfer_pma2/dm2zp.m b/+tfer_pma/dm2zp.m similarity index 100% rename from +tfer_pma2/dm2zp.m rename to +tfer_pma/dm2zp.m diff --git a/+tfer_pma2/get_setpoint.m b/+tfer_pma/get_setpoint.m similarity index 100% rename from +tfer_pma2/get_setpoint.m rename to +tfer_pma/get_setpoint.m diff --git a/+tfer_pma2/mp2zp.m b/+tfer_pma/mp2zp.m similarity index 100% rename from +tfer_pma2/mp2zp.m rename to +tfer_pma/mp2zp.m diff --git a/+tfer_pma2/prop_PMA.m b/+tfer_pma/prop_PMA.m similarity index 100% rename from +tfer_pma2/prop_PMA.m rename to +tfer_pma/prop_PMA.m diff --git a/+tfer_pma2/tfer_1C.m b/+tfer_pma/tfer_1C.m similarity index 100% rename from +tfer_pma2/tfer_1C.m rename to +tfer_pma/tfer_1C.m diff --git a/+tfer_pma2/tfer_1C_diff.m b/+tfer_pma/tfer_1C_diff.m similarity index 100% rename from +tfer_pma2/tfer_1C_diff.m rename to +tfer_pma/tfer_1C_diff.m diff --git a/+tfer_pma2/tfer_1C_pb.m b/+tfer_pma/tfer_1C_pb.m similarity index 100% rename from +tfer_pma2/tfer_1C_pb.m rename to +tfer_pma/tfer_1C_pb.m diff --git a/+tfer_pma2/tfer_1S.m b/+tfer_pma/tfer_1S.m similarity index 100% rename from +tfer_pma2/tfer_1S.m rename to +tfer_pma/tfer_1S.m diff --git a/+tfer_pma2/tfer_1S_diff.m b/+tfer_pma/tfer_1S_diff.m similarity index 100% rename from +tfer_pma2/tfer_1S_diff.m rename to +tfer_pma/tfer_1S_diff.m diff --git a/+tfer_pma2/tfer_1S_pb.m b/+tfer_pma/tfer_1S_pb.m similarity index 100% rename from +tfer_pma2/tfer_1S_pb.m rename to +tfer_pma/tfer_1S_pb.m diff --git a/+tfer_pma2/tfer_2C.m b/+tfer_pma/tfer_2C.m similarity index 100% rename from +tfer_pma2/tfer_2C.m rename to +tfer_pma/tfer_2C.m diff --git a/+tfer_pma2/tfer_2C_diff.m b/+tfer_pma/tfer_2C_diff.m similarity index 100% rename from +tfer_pma2/tfer_2C_diff.m rename to +tfer_pma/tfer_2C_diff.m diff --git a/+tfer_pma2/tfer_2S.m b/+tfer_pma/tfer_2S.m similarity index 100% rename from +tfer_pma2/tfer_2S.m rename to +tfer_pma/tfer_2S.m diff --git a/+tfer_pma2/tfer_2S_diff.m b/+tfer_pma/tfer_2S_diff.m similarity index 100% rename from +tfer_pma2/tfer_2S_diff.m rename to +tfer_pma/tfer_2S_diff.m diff --git a/+tfer_pma2/tfer_Ehara.m b/+tfer_pma/tfer_Ehara.m similarity index 100% rename from +tfer_pma2/tfer_Ehara.m rename to +tfer_pma/tfer_Ehara.m diff --git a/+tfer_pma2/tfer_FD.m b/+tfer_pma/tfer_FD.m similarity index 100% rename from +tfer_pma2/tfer_FD.m rename to +tfer_pma/tfer_FD.m diff --git a/+tfer_pma2/tfer_GE.m b/+tfer_pma/tfer_GE.m similarity index 100% rename from +tfer_pma2/tfer_GE.m rename to +tfer_pma/tfer_GE.m diff --git a/+tfer_pma2/tfer_GE_diff.m b/+tfer_pma/tfer_GE_diff.m similarity index 100% rename from +tfer_pma2/tfer_GE_diff.m rename to +tfer_pma/tfer_GE_diff.m diff --git a/+tfer_pma2/tfer_W1.m b/+tfer_pma/tfer_W1.m similarity index 100% rename from +tfer_pma2/tfer_W1.m rename to +tfer_pma/tfer_W1.m diff --git a/+tfer_pma2/tfer_W1_diff.m b/+tfer_pma/tfer_W1_diff.m similarity index 100% rename from +tfer_pma2/tfer_W1_diff.m rename to +tfer_pma/tfer_W1_diff.m diff --git a/+tfer_pma2/tfer_W1_pb.m b/+tfer_pma/tfer_W1_pb.m similarity index 100% rename from +tfer_pma2/tfer_W1_pb.m rename to +tfer_pma/tfer_W1_pb.m diff --git a/+tfer_pma2/tfer_tri.m b/+tfer_pma/tfer_tri.m similarity index 100% rename from +tfer_pma2/tfer_tri.m rename to +tfer_pma/tfer_tri.m diff --git a/+tfer_pma2/tridiag.m b/+tfer_pma/tridiag.m similarity index 100% rename from +tfer_pma2/tridiag.m rename to +tfer_pma/tridiag.m diff --git a/main_ehara.m b/main_ehara.m deleted file mode 100644 index 0482bb0..0000000 --- a/main_ehara.m +++ /dev/null @@ -1,81 +0,0 @@ - -% MAIN Script used in plotting different transfer functions. -% Author: Timothy Sipkens, 2019-06-25 -%=========================================================================% - - -%-- Initialize script ----------------------------------------------------% -clear; -close all; - -V = 10; % voltage to replicate Ehara et al. -omega = 1000*0.1047; % angular speed, converted from rpm to rad/s - -e = 1.60218e-19; % electron charge [C] -m = linspace(1e-10,7,601)*e; % vector of mass -% m = linspace(1e-10,6,601)*e; % vector of mass - -z = 1; % integer charge state - -rho_eff = 900; % effective density -d = (6.*m./(rho_eff.*pi)).^(1/3); - % specify mobility diameter vector with constant effective density - -prop = tfer_pma.prop_PMA('Ehara'); % get properties of the CPMA - -%=========================================================================% -%-- Finite difference solution -------------------------------------------% -tic; -[tfer_FD,~,n] = tfer_pma.tfer_FD([],... - m,d,1,prop,'V',V,'omega',omega); -t(1) = toc; - - -%=========================================================================% -%-- Transfer functions for different cases -------------------------------% -%-- Setup for centriputal force ------------------------------------------% -B = tfer_pma.dm2zp(d,z,prop.T,prop.p); -tau = B.*m; -D = prop.D(B); -sig = sqrt(2.*prop.L.*D./prop.v_bar); -D0 = D.*prop.L/(prop.del^2*prop.v_bar); % dimensionless diffusion coeff. - - -%-- Particle tracking approaches -----------------------------------------% -%-- Plug flow ------------------------------------------------------------% -%-- Method 1S ------------------------------% -tic; -[tfer_1S,G0_1S] = tfer_pma.tfer_1S([],m,d,z,prop,'V',V,'omega',omega); -t(2) = toc; - -%-- Method 1S, Ehara et al. ----------------% -tfer_Ehara = tfer_pma.tfer_Ehara([],m,d,z,prop,'V',V,'omega',omega); - - - -%-- Parabolic flow -------------------------------------------------------% -%-- Method 1S ------------------------------% -tic; -[tfer_1S_pb,G0_1S_pb] = tfer_pma.tfer_1S_pb([],m,d,z,prop,'V',V,'omega',omega); -t(8) = toc; - - - -%=========================================================================% -%-- Plot different transfer functions with respect to m/m* ---------------% -m_plot = m./e; - -figure(2); -plot(m_plot,tfer_1S); -hold on; -plot(m_plot,tfer_Ehara); -plot(m_plot,tfer_1S_pb); -plot(m_plot,min(tfer_FD,1),'k'); -hold off; - -% ylim([0,1.2]); - -xlabel('s') -ylabel('{\Lambda}') - - diff --git a/main_Ehara.m b/main_ehara2.m similarity index 100% rename from main_Ehara.m rename to main_ehara2.m From 6c118e4cd638a5be2fc9d2a4f07db6003c12e059 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 19 Sep 2019 15:58:05 -0700 Subject: [PATCH 20/22] Revert back to original naming (pt 2) --- main_ehara2.m => main_ehara.m | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename main_ehara2.m => main_ehara.m (100%) diff --git a/main_ehara2.m b/main_ehara.m similarity index 100% rename from main_ehara2.m rename to main_ehara.m From feba447b722ce9932615bd93d304e2f12adbd943 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Thu, 19 Sep 2019 16:01:26 -0700 Subject: [PATCH 21/22] Update README.md --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index d7228c6..68a86c4 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# mat-tfer-PMA +# mat-tfer-pma The attached MATLAB functions and scripts are intended to reproduce the results of the associated paper (submitted). They evaluate the transfer @@ -37,7 +37,7 @@ The functions share common inputs: evaluated), 5. *prop* - a struct that contains the properties of the particle mass analyzer - (a sample script to generate this quantity is include as `prop_CPMA.m`), and + (a sample script to generate this quantity is include as `prop_PMA.m`), and 6. *varargin* (optional) - name-value pairs to specify either the equivalent resolution, inner electrode angular speed, or voltage. From 551c92c8317d00c0709524a0e225a0312c9bb735 Mon Sep 17 00:00:00 2001 From: "Timothy A. Sipkens" Date: Fri, 20 Sep 2019 16:49:23 -0700 Subject: [PATCH 22/22] Updates to get_setpoint.m GET_SETPOINT is now a function, cleaning warning in existing files. --- +tfer_pma/get_setpoint.m | 15 ++++++--------- +tfer_pma/tfer_1C.m | 5 +++-- +tfer_pma/tfer_1C_pb.m | 4 +++- +tfer_pma/tfer_1S.m | 4 +++- +tfer_pma/tfer_1S_pb.m | 4 +++- +tfer_pma/tfer_2C.m | 4 +++- +tfer_pma/tfer_2S.m | 4 +++- +tfer_pma/tfer_Ehara.m | 4 +++- +tfer_pma/tfer_FD.m | 5 +++-- +tfer_pma/tfer_GE.m | 4 +++- +tfer_pma/tfer_W1.m | 4 +++- +tfer_pma/tfer_W1_pb.m | 4 +++- +tfer_pma/tfer_tri.m | 3 ++- 13 files changed, 41 insertions(+), 23 deletions(-) diff --git a/+tfer_pma/get_setpoint.m b/+tfer_pma/get_setpoint.m index 80e6044..46ccd47 100644 --- a/+tfer_pma/get_setpoint.m +++ b/+tfer_pma/get_setpoint.m @@ -2,7 +2,8 @@ % GET_SETPOINT Script to evaluate setpoint parameters including C0, alpha, and beta. % Author: Timothy A. Sipkens, 2019-05-01 %=========================================================================% -% + +function [sp,tau,C0,D] = get_setpoint(m_star,m,d,z,prop,varargin) %-------------------------------------------------------------------------% % Required variables: % m_star Mass corresponding to the measurement set point of the APM @@ -28,14 +29,10 @@ %-- Parse inputs ---------------------------------------------------------% -if ~exist('z','var'); z = []; end if isempty(z); z = 1; end % if integer charge is not specified, use z = 1 -if ~exist('m_star','var'); m_star = []; end - - %-- Parse inputs for setpoint --------------------------------------------% -if exist('varargin','var') % parse input name-value pair, if it exists +if ~isempty(varargin) % parse input name-value pair, if it exists if length(varargin)==2 sp.(varargin{1}) = varargin{2}; elseif length(varargin)==4 @@ -53,7 +50,7 @@ e = 1.60218e-19; % electron charge [C] q = z.*e; % particle charge -if ~exist('d','var') % evaluate mechanical mobility +if isempty(d) % evaluate mechanical mobility warning('Invoking mass-mobility relation to determine Zp.'); B = tfer_pma.mp2zp(m,z,prop.T,prop.p); else @@ -157,7 +154,7 @@ % copy center mass to sp % center mass is for a singly charged particle -B_star = tfer_pma.mp2zp(m_star,1,prop.T,prop.p); +% B_star = tfer_pma.mp2zp(m_star,1,prop.T,prop.p); C0 = sp.V.*q./log(1/prop.r_hat); % calcualte recurring C0 parameter - +end diff --git a/+tfer_pma/tfer_1C.m b/+tfer_pma/tfer_1C.m index ec4ab2b..4b39f23 100644 --- a/+tfer_pma/tfer_1C.m +++ b/+tfer_pma/tfer_1C.m @@ -21,8 +21,9 @@ % G0 Function mapping final to initial radial position %-------------------------------------------------------------------------% - -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- 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)); diff --git a/+tfer_pma/tfer_1C_pb.m b/+tfer_pma/tfer_1C_pb.m index 233c30b..a5f9563 100644 --- a/+tfer_pma/tfer_1C_pb.m +++ b/+tfer_pma/tfer_1C_pb.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- 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)); diff --git a/+tfer_pma/tfer_1S.m b/+tfer_pma/tfer_1S.m index caaa97e..2666dd1 100644 --- a/+tfer_pma/tfer_1S.m +++ b/+tfer_pma/tfer_1S.m @@ -21,7 +21,9 @@ % G0 Function mapping final to initial radial position %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_1S_pb.m b/+tfer_pma/tfer_1S_pb.m index d0eefdf..a80563f 100644 --- a/+tfer_pma/tfer_1S_pb.m +++ b/+tfer_pma/tfer_1S_pb.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_2C.m b/+tfer_pma/tfer_2C.m index 0b75275..696d7da 100644 --- a/+tfer_pma/tfer_2C.m +++ b/+tfer_pma/tfer_2C.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- 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)); diff --git a/+tfer_pma/tfer_2S.m b/+tfer_pma/tfer_2S.m index 06ecc44..bfbb879 100644 --- a/+tfer_pma/tfer_2S.m +++ b/+tfer_pma/tfer_2S.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_Ehara.m b/+tfer_pma/tfer_Ehara.m index f880d71..ee9f1a2 100644 --- a/+tfer_pma/tfer_Ehara.m +++ b/+tfer_pma/tfer_Ehara.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_FD.m b/+tfer_pma/tfer_FD.m index aa077d1..ff5c284 100644 --- a/+tfer_pma/tfer_FD.m +++ b/+tfer_pma/tfer_FD.m @@ -27,8 +27,9 @@ % laboratory. %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) - +[sp,tau,C0,D] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Discretize the space -------------------------------------------------% nr = 200; diff --git a/+tfer_pma/tfer_GE.m b/+tfer_pma/tfer_GE.m index 9d0b3fe..184d57d 100644 --- a/+tfer_pma/tfer_GE.m +++ b/+tfer_pma/tfer_GE.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_W1.m b/+tfer_pma/tfer_W1.m index b9751a7..d3d0c6e 100644 --- a/+tfer_pma/tfer_W1.m +++ b/+tfer_pma/tfer_W1.m @@ -22,7 +22,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_W1_pb.m b/+tfer_pma/tfer_W1_pb.m index b84ecc6..262ebf9 100644 --- a/+tfer_pma/tfer_W1_pb.m +++ b/+tfer_pma/tfer_W1_pb.m @@ -23,7 +23,9 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint (parses d and z) +[sp,tau,C0] = ... + tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) %-- Estimate equilibrium radius ------------------------------------------% if round((sqrt(C0./m_star)-sqrt(C0./m_star-4*sp.alpha*sp.beta))/(2*sp.alpha),15)==prop.rc diff --git a/+tfer_pma/tfer_tri.m b/+tfer_pma/tfer_tri.m index 69ca2da..a732b4c 100644 --- a/+tfer_pma/tfer_tri.m +++ b/+tfer_pma/tfer_tri.m @@ -22,7 +22,8 @@ %-------------------------------------------------------------------------% -tfer_pma.get_setpoint; % get setpoint +sp = tfer_pma.get_setpoint(m_star,m,d,z,prop,varargin{:}); + % get setpoint (parses d and z) if ~isfield(sp,'m_max') % if m_max was not specified n_B = -0.6436;