Merge pull request #7 from tsipkens/developer

Developer

+tfer_PMA/G_fun.m → +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) <r1
+        if rL(jj)<r1
         % particles cannot exist here, return r1
             G(jj,ii) = r1;
 
@@ -42,29 +43,27 @@
             G(jj,ii) = r2;
 
         elseif and(rL(jj)<(rs(ii)+5*offset),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)<rs(ii) % should be equivalent to else
         % if rL < rs, then r0 < rL
 
-            if condit(ii)
+            if condit(ii) % zero is in interval [r1,rL]
                 if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),r1,ii))
                 % if sign does not change in interval [r1,rL], return r1
                     G(jj,ii) = r1;
 
-                else
-                % find zero in interval [r1,rL]
+                else % find zero in interval [r1,rL]
                     G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [r1,rL(jj)]);
 
                 end
 
-            else
+            else % zero is in interval [rL,min(r2,rs)]
                 if sign(min_fun(rL(jj),rL(jj),ii))==sign(min_fun(rL(jj),min(rs(ii)-offset,r2),ii))
                 % if sign does not change in interval [rL,min(r2,rs)], return min(r2,rs)
                     G(jj,ii) = min(rs(ii)-offset,r2);
 
-                else
-                % find zero in interval [rL,min(r2,rs)]
+                else % find zero in interval [rL,min(r2,rs)]
                     G(jj,ii) = fzero(@(r0) min_fun(rL(jj),r0,ii), [rL(jj),min(rs(ii)-offset,r2)]);
 
                 end

+tfer_PMA/dm2zp.m → +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

+tfer_PMA/get_setpoint.m → +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
@@ -15,7 +16,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.)
@@ -27,27 +28,16 @@
 %-------------------------------------------------------------------------%
 
 
-%-- Set up mobility calculations and parse inputs ------------------------%
-if ~exist('z','var') % if integer charge is not specified, use z = 1
-    z = 1;
-end
-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 ---------------------------------------------------------%
+if isempty(z); z = 1; end % if integer charge is not specified, use z = 1
 
 %-- 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
+        sp.(varargin{1}) = varargin{2};
+        sp.(varargin{3}) = varargin{4};
     else
         sp.Rm = 3; % by default use resolution, with value of 3
     end
@@ -56,33 +46,90 @@
 end
 
 
-%-- Proceed depnding on which setpoint parameter is specified ------------%
-if isfield(sp,'omega1') % if angular speed of inner electrode is specified
+%-- Set up mobility calculations -----------------------------------------%
+e = 1.60218e-19; % electron charge [C]
+q = z.*e; % particle charge
+
+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
+    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.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);
+
+    %-- 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 ------------------------%
     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
 
+    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)));
+        % 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;
 
     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
 
     %-- 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)));
@@ -95,15 +142,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
 
-
+end

+tfer_PMA/mp2zp.m → +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

+tfer_PMA/prop_CPMA.m → +tfer_pma/prop_PMA.m

@@ -1,9 +1,9 @@
 
 % 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)
+function [prop] = prop_PMA(opt)
 %-------------------------------------------------------------------------%
 % Input:
 %   opt         Options string specifying parameter set
@@ -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')
@@ -35,13 +35,42 @@
     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 = 0.5/1000/60; % aerosol flowrate [m^3/s], assumed
+    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]
+
+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