-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathget_setpoint.m
247 lines (172 loc) · 9.12 KB
/
get_setpoint.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
% GET_SETPOINT Generate setpoint parameter structure from available parameters.
% Author: Timothy Sipkens, 2019-11-21
%
% Required variables:
% prop Properties of particle mass analyzer
% varargin Name-value pairs for setpoint
% (two values required, if one value is specified, uses Rm = 3)
% ('m_star',double) - Setpoint mass in kg (i.e., 1e-18 for a 1 fg particle)
% ('Rm',double) - Resolution
% ('omega1',double) - Angular speed of inner electrode
% ('V',double) - Setpoint voltage
%
% Sample outputs:
% sp Struct containing mutliple setpoint parameters (V, alpha, etc.)
% m_star Setpoint mass, assuming a singly charged particle, in kg
%
% Notes:
% As a script, this code uses variables currently in the
% workspace. This script is also used to parse some of the inputs
% to the various transfer functions, including the existence of
% the integer charge state and particle mobility.
%=========================================================================%
%== GET_SETPOINT =========================================================%
% Wrapper function for get_setpoint0 to loop through a range of setpoints.
function [sp,m_star] = get_setpoint(prop,varargin)
%-- Parse inputs ---------------------------------------------------------%
if ~exist('prop','var'); prop = struct(); end
if isempty(prop); prop = prop_pma; end
if isempty(varargin) % if empty input, return empty structure
sp.m_star = []; sp.V = []; sp.Rm = [];
sp.omega = []; sp.omega1 = []; sp.omega2 = [];
sp.alpha = []; sp.beta = [];
sp.m_max = [];
m_star = [];
return;
end
if length(varargin)==2; varargin = [varargin,'Rm',3]; end % by default use Rm = 3
%-------------------------------------------------------------------------%
n = max(length(varargin{2}),length(varargin{4})); % number of setpoints
if length(varargin{2})~=n; varargin{2} = varargin{2}.*ones(n,1); end
if length(varargin{4})~=n; varargin{4} = varargin{4}.*ones(n,1); end
% expand scalar entries to n setpoints
for ii=n:-1:1 % loop through prescribed setpoints
[sp(ii),m_star(ii)] = get_setpoint0(prop,...
varargin{1},varargin{2}(ii),...
varargin{3},varargin{4}(ii));
% get CPMA setpoint parameters for the iith setpoint
end
end
%== GET_SETPOINT0 ========================================================%
% Get parameters for a single setpoint (original function).
function [sp,m_star] = get_setpoint0(prop,varargin)
%-- Initialize sp structure ----------------------------------------------%
sp = struct('m_star',[],'V',[],'Rm',[],'omega',[],...
'omega1',[],'omega2',[],'alpha',[],'beta',[],'m_max',[]);
% default empty structure
% sp.m_star corresponds the mass at the given setpoint
% for a singly-charged particle
%-- Parse inputs ---------------------------------------------------------%
for ii=1:2:length(varargin)
sp.(varargin{ii}) = varargin{ii+1};
end
%-------------------------------------------------------------------------%
e = 1.60218e-19; % electron charge [C]
%== Proceed depending on which setpoint parameters are specified =========%
%== CASE 1: m_star is not specified, use V and omega =====================%
if isempty(sp.m_star)
% case if m_star is not specified
% requires that voltage, 'V', and speed, 'omega' are specified
%-- Evaluate angular speed of inner electrode ------------------------%
if isempty(sp.omega1)
sp.omega1 = sp.omega./...
((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);
end
%-- 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.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.omega = sp.alpha+sp.beta./(prop.rc.^2);
sp.omega2 = sp.alpha+sp.beta./(prop.r2.^2);
%== CASE 2: m_star and omega1 are specified ==============================%
elseif ~isempty(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 = sp.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);
%== CASE 3: m_star and omega are specified ===============================%
elseif ~isempty(sp.omega) % if angular speed at gap center is specified
%-- Evaluate angular speed of inner electrode ------------------------%
sp.omega1 = sp.omega./...
((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);
%-- 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 = sp.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);
%== CASE 4: m_star and V are specified ===================================%
elseif ~isempty(sp.V) % if voltage is specified
v_theta_rc = sqrt(sp.V.*e./(sp.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);
%== CASE 5: m_star and Rm are specified ==================================%
elseif ~isempty(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 = get_nb(sp.m_star,prop);
B_star = mp2zp(sp.m_star,1,prop.T,prop.p,prop);
% involves invoking mass-mobility relation
% z = 1 for the setpoint
sp.m_max = sp.m_star*(1/sp.Rm+1);
sp.omega = sqrt(prop.Q/(sp.m_star*B_star*2*pi*prop.rc^2*prop.L*...
((sp.m_max/sp.m_star)^(n_B+1)-(sp.m_max/sp.m_star)^n_B)));
%-- Evaluate angular speed of inner electrode ------------------------%
sp.omega1 = sp.omega./...
((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);
%-- 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.V = sp.m_star.*log(1/prop.r_hat)./e.*(sp.alpha.*prop.rc+sp.beta./prop.rc).^2;
else % other cases are not supported, output error
error('Invalid setpoint parameters specified.');
end
%-- Calculate resolution -------------------------------------------------%
if isempty(sp.Rm) % if resolution is not specified
[sp.Rm,sp.m_max] = get_resolution(sp.m_star,sp.omega,prop);
% evaluate resolution in corresponding subfunction
% involves a minimization routine
end
m_star = sp.m_star; % output m_star independently
end
%== GET_RESOLUTION =======================================================%
% Solver to evaluate the resolution from m_star and prop.
function [Rm,m_max] = get_resolution(m_star,omega,prop)
n_B = get_nb(m_star,prop);
B_star = mp2zp(m_star,1,...
prop.T,prop.p,prop); % mechanical mobility for z = 1
t0 = prop.Q/(m_star*B_star*2*pi*prop.L*...
omega^2*prop.rc^2); % RHS of Eq. (10) in Reveall et al.
m_rat = @(Rm) 1/Rm+1; % function for mass ratio
fun = @(Rm) (m_rat(Rm))^(n_B+1)-(m_rat(Rm))^n_B; % LHS of Eq. (10) in Reveall et al.
Rm = fminsearch(@(Rm) (t0-fun(Rm))^2,5); % minimization ot find Rm
m_max = m_star*(1/Rm+1); % approx. upper end of non-diffusing tfer. function
end
%== GET_NB ===============================================================%
% Function to evaluate n_B constant. Taken from Olfert laboratory.
% Note: Previous versions of this program would output a constant
% value of n_B = -0.6436. This will cause some rather minor
% compatiblity issues.
function n_B = get_nb(m_star,prop)
m_high = m_star*1.001; % perturb m_star up
m_low = m_star*.999; % perturb m_star down
B_high = mp2zp(m_high,1,prop.T,prop.p,prop);
B_low = mp2zp(m_low,1,prop.T,prop.p,prop);
n_B = log10(B_high/B_low)/log10(m_high/m_low); % constant from Reveall et al.
% n_B = -0.6436; % deprecated value
end