1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion = 6.0 |
---|
3 | |
---|
4 | //CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
5 | //CCCCCCCC |
---|
6 | //C SUBROUTINE FOR THE CALCULATION OF THE SCATTERING FUNCTIONS |
---|
7 | //C OF RODLIKE MICELLES. METHODLOGY FOLLOWS THAT OF PEDERSEN AND |
---|
8 | //C SCHURTENBERGER, MACORMOLECULES, VOL 29,PG 7602, 1996. |
---|
9 | //C WITH EXCULDED VOLUME EFFECTS (METHOD 3) |
---|
10 | // |
---|
11 | // - copied directly from FORTRAN code supplied by Jan Pedersen |
---|
12 | // SRK - 2002, but shows discontinuity at Qlb = 3.1 |
---|
13 | // |
---|
14 | // Jan 2006 - re-worked FORTRAN correcting typos in paper: now is smooth, but |
---|
15 | // the splicing is actually at Qlb = 2, which is not what the paper |
---|
16 | // says is to be done (and not from earlier models) |
---|
17 | // |
---|
18 | // July 2006 - now is CORRECT with Wei-Ren's changes to the code |
---|
19 | // Matlab code was not too difficult to convert to Igor (only a few hours...) |
---|
20 | // |
---|
21 | Proc PlotFlexExclVolCyl(num,qmin,qmax) |
---|
22 | Variable num=128,qmin=0.001,qmax=0.7 |
---|
23 | Prompt num "Enter number of data points for model: " |
---|
24 | Prompt qmin "Enter minimum q-value (^-1) for model: " |
---|
25 | Prompt qmax "Enter maximum q-value (^-1) for model: " |
---|
26 | |
---|
27 | Make/O/D/n=(num) xwave_fle,ywave_fle |
---|
28 | xwave_fle = alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) |
---|
29 | Make/O/D coef_fle = {1.,1000,100,20,1.0e-6,6.3e-6,0.0001} |
---|
30 | make/o/t parameters_fle = {"scale","Contour Length (A)","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"} |
---|
31 | Edit parameters_fle,coef_fle |
---|
32 | |
---|
33 | Variable/G root:g_fle |
---|
34 | g_fle := FlexExclVolCyl(coef_fle,ywave_fle,xwave_fle) |
---|
35 | Display ywave_fle vs xwave_fle |
---|
36 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
37 | Label bottom "q (\\S-1\\M)" |
---|
38 | Label left "Intensity (cm\\S-1\\M)" |
---|
39 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
40 | |
---|
41 | AddModelToStrings("FlexExclVolCyl","coef_fle","fle") |
---|
42 | End |
---|
43 | /////////////////////////////////////////////////////////// |
---|
44 | |
---|
45 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
---|
46 | Proc PlotSmearedFlexExclVolCyl(str) |
---|
47 | String str |
---|
48 | Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) |
---|
49 | |
---|
50 | // if any of the resolution waves are missing => abort |
---|
51 | if(ResolutionWavesMissingDF(str)) //updated to NOT use global strings (in GaussUtils) |
---|
52 | Abort |
---|
53 | endif |
---|
54 | |
---|
55 | SetDataFolder $("root:"+str) |
---|
56 | |
---|
57 | // Setup parameter table for model function |
---|
58 | Make/O/D smear_coef_fle = {1.,1000,100,20,1.0e-6,6.3e-6,0.0001} |
---|
59 | make/o/t smear_parameters_fle = {"scale","Contour Length (A)","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"} |
---|
60 | Edit smear_parameters_fle,smear_coef_fle |
---|
61 | |
---|
62 | // output smeared intensity wave, dimensions are identical to experimental QSIG values |
---|
63 | // make extra copy of experimental q-values for easy plotting |
---|
64 | Duplicate/O $(str+"_q") smeared_fle,smeared_qvals |
---|
65 | SetScale d,0,0,"1/cm",smeared_fle |
---|
66 | |
---|
67 | Variable/G gs_fle=0 |
---|
68 | gs_fle := fSmearedFlexExclVolCyl(smear_coef_fle,smeared_fle,smeared_qvals) //this wrapper fills the STRUCT |
---|
69 | |
---|
70 | Display smeared_fle vs smeared_qvals |
---|
71 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
72 | Label bottom "q (\\S-1\\M)" |
---|
73 | Label left "Intensity (cm\\S-1\\M)" |
---|
74 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
75 | |
---|
76 | SetDataFolder root: |
---|
77 | AddModelToStrings("SmearedFlexExclVolCyl","smear_coef_fle","fle") |
---|
78 | End |
---|
79 | |
---|
80 | |
---|
81 | //AAO version, uses XOP if available |
---|
82 | // simply calls the original single point calculation with |
---|
83 | // a wave assignment (this will behave nicely if given point ranges) |
---|
84 | Function FlexExclVolCyl(cw,yw,xw) : FitFunc |
---|
85 | Wave cw,yw,xw |
---|
86 | |
---|
87 | #if exists("FlexExclVolCylX") |
---|
88 | yw = FlexExclVolCylX(cw,xw) |
---|
89 | #else |
---|
90 | yw = fFlexExclVolCyl(cw,xw) |
---|
91 | #endif |
---|
92 | return(0) |
---|
93 | End |
---|
94 | |
---|
95 | // |
---|
96 | Function fFlexExclVolCyl(ww,x) |
---|
97 | Wave ww |
---|
98 | Variable x |
---|
99 | |
---|
100 | //nice names to the input params |
---|
101 | //ww[0] = scale |
---|
102 | //ww[1] = L [A] |
---|
103 | //ww[2] = B [A] |
---|
104 | //ww[3] = rad [A] cross-sectional radius |
---|
105 | //ww[4] = sld cylinder |
---|
106 | //ww[5] = sld solvent [A^-2] |
---|
107 | //ww[6] = bkg [cm-1] |
---|
108 | Variable scale,L,B,bkg,rad,qr,cont,sldc,slds |
---|
109 | |
---|
110 | scale = ww[0] |
---|
111 | L = ww[1] |
---|
112 | B = ww[2] |
---|
113 | rad = ww[3] |
---|
114 | sldc = ww[4] |
---|
115 | slds = ww[5] |
---|
116 | bkg = ww[6] |
---|
117 | qr = x*rad //used for cross section contribution only |
---|
118 | |
---|
119 | cont = sldc-slds |
---|
120 | |
---|
121 | Variable flex,crossSect |
---|
122 | flex = Sk_WR(x,L,B) |
---|
123 | |
---|
124 | crossSect = (2*bessJ(1,qr)/qr)^2 |
---|
125 | |
---|
126 | //normalize form factor by multiplying by cylinder volume * cont^2 |
---|
127 | // then convert to cm-1 by multiplying by 10^8 |
---|
128 | // then scale = phi |
---|
129 | |
---|
130 | flex *= crossSect |
---|
131 | flex *= Pi*rad*rad*L |
---|
132 | flex *= cont^2 |
---|
133 | flex *= 1.0e8 |
---|
134 | |
---|
135 | return (scale*flex + bkg) |
---|
136 | End |
---|
137 | |
---|
138 | //////////////////WRC corrected code below |
---|
139 | // main function |
---|
140 | function Sk_WR(q,L,b) |
---|
141 | Variable q,L,b |
---|
142 | // |
---|
143 | Variable p1,p2,p1short,p2short,q0,qconnect |
---|
144 | Variable c,epsilon,ans,q0short,Sexvmodify |
---|
145 | |
---|
146 | p1 = 4.12 |
---|
147 | p2 = 4.42 |
---|
148 | p1short = 5.36 |
---|
149 | p2short = 5.62 |
---|
150 | q0 = 3.1 |
---|
151 | qconnect = q0/b |
---|
152 | // |
---|
153 | q0short = max(1.9/sqrt(Rgsquareshort(q,L,b)),3) |
---|
154 | |
---|
155 | // |
---|
156 | if(L/b > 10) |
---|
157 | C = 3.06/(L/b)^0.44 |
---|
158 | epsilon = 0.176 |
---|
159 | else |
---|
160 | C = 1 |
---|
161 | epsilon = 0.170 |
---|
162 | endif |
---|
163 | // |
---|
164 | |
---|
165 | if( L > 4*b ) // Longer Chains |
---|
166 | if (q*b <= 3.1) |
---|
167 | //Modified by Yun on Oct. 15, |
---|
168 | Sexvmodify = Sexvnew(q, L, b) |
---|
169 | ans = Sexvmodify + C * (4/15 + 7./(15*u_WR(q,L,b)) - (11/15 + 7./(15*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L) |
---|
170 | else //q(i)*b > 3.1 |
---|
171 | ans = a1long(q, L, b, p1, p2, q0)/((q*b)^p1) + a2long(q, L, b, p1, p2, q0)/((q*b)^p2) + pi/(q*L) |
---|
172 | endif |
---|
173 | else //L <= 4*b Shorter Chains |
---|
174 | if (q*b <= max(1.9/sqrt(Rgsquareshort(q,L,b)),3) ) |
---|
175 | if (q*b<=0.01) |
---|
176 | ans = 1 - Rgsquareshort(q,L,b)*(q^2)/3 |
---|
177 | else |
---|
178 | ans = Sdebye1(q,L,b) |
---|
179 | endif |
---|
180 | else //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) |
---|
181 | ans = a1short(q,L,b,p1short,p2short,q0short)/((q*b)^p1short) + a2short(q,L,b,p1short,p2short,q0short)/((q*b)^p2short) + pi/(q*L) |
---|
182 | endif |
---|
183 | endif |
---|
184 | |
---|
185 | return(ans) |
---|
186 | end |
---|
187 | |
---|
188 | //WR named this w (too generic) |
---|
189 | Function w_WR(x) |
---|
190 | Variable x |
---|
191 | |
---|
192 | //C4 = 1.523; |
---|
193 | //C5 = 0.1477; |
---|
194 | Variable yy |
---|
195 | yy = 0.5*(1 + tanh((x - 1.523)/0.1477)) |
---|
196 | |
---|
197 | return (yy) |
---|
198 | end |
---|
199 | |
---|
200 | // |
---|
201 | function u1(q,L,b) |
---|
202 | Variable q,L,b |
---|
203 | Variable yy |
---|
204 | |
---|
205 | yy = Rgsquareshort(q,L,b)*q^2 |
---|
206 | |
---|
207 | return yy |
---|
208 | end |
---|
209 | |
---|
210 | // was named u |
---|
211 | function u_WR(q,L,b) |
---|
212 | Variable q,L,b |
---|
213 | |
---|
214 | Variable yy |
---|
215 | yy = Rgsquare(q,L,b)*(q^2) |
---|
216 | return yy |
---|
217 | end |
---|
218 | |
---|
219 | |
---|
220 | |
---|
221 | // |
---|
222 | function Rgsquarezero(q,L,b) |
---|
223 | Variable q,L,b |
---|
224 | |
---|
225 | Variable yy |
---|
226 | yy = (L*b/6) * (1 - 1.5*(b/L) + 1.5*(b/L)^2 - 0.75*(b/L)^3*(1 - exp(-2*(L/b)))) |
---|
227 | |
---|
228 | return yy |
---|
229 | end |
---|
230 | |
---|
231 | // |
---|
232 | function Rgsquareshort(q,L,b) |
---|
233 | Variable q,L,b |
---|
234 | |
---|
235 | Variable yy |
---|
236 | yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b) |
---|
237 | |
---|
238 | return yy |
---|
239 | end |
---|
240 | |
---|
241 | // |
---|
242 | function Rgsquare(q,L,b) |
---|
243 | Variable q,L,b |
---|
244 | |
---|
245 | Variable yy |
---|
246 | yy = AlphaSquare(L/b)*L*b/6 |
---|
247 | |
---|
248 | return yy |
---|
249 | end |
---|
250 | |
---|
251 | // |
---|
252 | function AlphaSquare(x) |
---|
253 | Variable x |
---|
254 | |
---|
255 | Variable yy |
---|
256 | yy = (1 + (x/3.12)^2 + (x/8.67)^3)^(0.176/3) |
---|
257 | |
---|
258 | return yy |
---|
259 | end |
---|
260 | |
---|
261 | // |
---|
262 | function miu(x) |
---|
263 | Variable x |
---|
264 | |
---|
265 | Variable yy |
---|
266 | yy = (1/8)*(9*x - 2 + 2*log(1 + x)/x)*exp(1/2.565*(1/x + (1 - 1/(x^2))*log(1 + x))) |
---|
267 | |
---|
268 | return yy |
---|
269 | end |
---|
270 | |
---|
271 | // |
---|
272 | function Sdebye(q,L,b) |
---|
273 | Variable q,L,b |
---|
274 | |
---|
275 | Variable yy |
---|
276 | yy = 2*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1)/((u_WR(q,L,b))^2) |
---|
277 | |
---|
278 | return yy |
---|
279 | end |
---|
280 | |
---|
281 | // |
---|
282 | function Sdebye1(q,L,b) |
---|
283 | Variable q,L,b |
---|
284 | |
---|
285 | Variable yy |
---|
286 | yy = 2*(exp(-u1(q,L,b)) + u1(q,L,b) -1)/((u1(q,L,b))^2) |
---|
287 | |
---|
288 | return yy |
---|
289 | end |
---|
290 | |
---|
291 | // |
---|
292 | function Sexv(q,L,b) |
---|
293 | Variable q,L,b |
---|
294 | |
---|
295 | Variable yy,C1,C2,C3,miu,Rg2 |
---|
296 | C1=1.22 |
---|
297 | C2=0.4288 |
---|
298 | C3=-1.651 |
---|
299 | miu = 0.585 |
---|
300 | |
---|
301 | Rg2 = Rgsquare(q,L,b) |
---|
302 | |
---|
303 | yy = (1 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) + C2*(q*sqrt(Rg2))^(-2/miu) + C3*(q*sqrt(Rg2))^(-3/miu)) |
---|
304 | |
---|
305 | return yy |
---|
306 | end |
---|
307 | |
---|
308 | // this must be WR modified version |
---|
309 | function Sexvnew(q,L,b) |
---|
310 | Variable q,L,b |
---|
311 | |
---|
312 | Variable yy,C1,C2,C3,miu |
---|
313 | C1=1.22 |
---|
314 | C2=0.4288 |
---|
315 | C3=-1.651 |
---|
316 | miu = 0.585 |
---|
317 | |
---|
318 | //calculating the derivative to decide on the corection (cutoff) term? |
---|
319 | // I have modified this from WRs original code |
---|
320 | Variable del=1.05,C_star2,Rg2 |
---|
321 | if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0 ) |
---|
322 | C_star2 = 0 |
---|
323 | else |
---|
324 | C_star2 = 1 |
---|
325 | endif |
---|
326 | |
---|
327 | Rg2 = Rgsquare(q,L,b) |
---|
328 | |
---|
329 | yy = (1 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) + C2*(q*sqrt(Rg2))^(-2/miu) + C3*(q*sqrt(Rg2))^(-3/miu)) |
---|
330 | |
---|
331 | return yy |
---|
332 | end |
---|
333 | |
---|
334 | |
---|
335 | |
---|
336 | // these are the messy ones |
---|
337 | function a2short(q, L, b, p1short, p2short, q0) |
---|
338 | Variable q, L, b, p1short, p2short, q0 |
---|
339 | |
---|
340 | Variable yy,Rg2_sh |
---|
341 | Rg2_sh = Rgsquareshort(q,L,b) |
---|
342 | |
---|
343 | Variable t1 |
---|
344 | t1 = ((q0^2*Rg2_sh)/b^2) |
---|
345 | |
---|
346 | //E is the number e |
---|
347 | yy = ((-(1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^(-4 + p2short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p1short + 2*b^3*E^t1*L*p1short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p1short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p1short*pi*q0^3*Rg2_sh^2))))))) |
---|
348 | |
---|
349 | return yy |
---|
350 | end |
---|
351 | |
---|
352 | // |
---|
353 | function a1short(q, L, b, p1short, p2short, q0) |
---|
354 | Variable q, L, b, p1short, p2short, q0 |
---|
355 | |
---|
356 | Variable yy,Rg2_sh |
---|
357 | Rg2_sh = Rgsquareshort(q,L,b) |
---|
358 | |
---|
359 | Variable t1 |
---|
360 | t1 = ((q0^2*Rg2_sh)/b^2) |
---|
361 | |
---|
362 | yy = ((1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^((-4) + p1short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p2short + 2*b^3*E^t1*L*p2short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p2short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p2short*pi*q0^3*Rg2_sh^2)))))) |
---|
363 | |
---|
364 | return yy |
---|
365 | end |
---|
366 | |
---|
367 | // this one will be lots of trouble |
---|
368 | function a2long(q, L, b, p1, p2, q0) |
---|
369 | variable q, L, b, p1, p2, q0 |
---|
370 | |
---|
371 | Variable yy,c1,c2,c3,c4,c5,miu,c,Rg2,rRg |
---|
372 | Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13 |
---|
373 | |
---|
374 | if( L/b > 10) |
---|
375 | C = 3.06/(L/b)^0.44 |
---|
376 | else |
---|
377 | C = 1 |
---|
378 | endif |
---|
379 | |
---|
380 | C1 = 1.22 |
---|
381 | C2 = 0.4288 |
---|
382 | C3 = -1.651 |
---|
383 | C4 = 1.523 |
---|
384 | C5 = 0.1477 |
---|
385 | miu = 0.585 |
---|
386 | |
---|
387 | Rg2 = Rgsquare(q,L,b) |
---|
388 | t1 = (1/(b* p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2))) |
---|
389 | t2 = (b*C*(((-1*((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b)))/L |
---|
390 | t3 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(2*C5) |
---|
391 | t4 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(C5*q0^4*Rg2^2) |
---|
392 | t5 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2) |
---|
393 | t6 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^5*Rg2^2) |
---|
394 | t7 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu)) |
---|
395 | t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
---|
396 | t9 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2))))/L |
---|
397 | t10 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2) |
---|
398 | |
---|
399 | |
---|
400 | yy = ((-1*(t1* (((-q0^(-p1))*(((b^2*pi)/(L*q0^2) + t2 + t3 - t4 + t5 - t6 + 1/2*t7*t8)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) + t9 + t10 + 1/2*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))) |
---|
401 | |
---|
402 | |
---|
403 | return yy |
---|
404 | end |
---|
405 | |
---|
406 | //need to define this on my own |
---|
407 | Function sech_WR(x) |
---|
408 | variable x |
---|
409 | |
---|
410 | return(1/cosh(x)) |
---|
411 | end |
---|
412 | // |
---|
413 | function a1long(q, L, b, p1, p2, q0) |
---|
414 | Variable q, L, b, p1, p2, q0 |
---|
415 | |
---|
416 | Variable yy,c,c1,c2,c3,c4,c5,miu,Rg2,rRg |
---|
417 | Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16 |
---|
418 | |
---|
419 | if( L/b > 10) |
---|
420 | C = 3.06/(L/b)^0.44 |
---|
421 | else |
---|
422 | C = 1 |
---|
423 | endif |
---|
424 | |
---|
425 | C1 = 1.22 |
---|
426 | C2 = 0.4288 |
---|
427 | C3 = -1.651 |
---|
428 | C4 = 1.523 |
---|
429 | C5 = 0.1477 |
---|
430 | miu = 0.585 |
---|
431 | |
---|
432 | Rg2 = Rgsquare(q,L,b) |
---|
433 | t1 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2)))) |
---|
434 | t2 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
---|
435 | t3 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu))) |
---|
436 | t4 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
---|
437 | t5 = (1/(b*p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2))) |
---|
438 | t6 = (b*C*(((-((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b))) |
---|
439 | t7 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2) |
---|
440 | t8 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2) |
---|
441 | t9 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
---|
442 | t10 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
---|
443 | t11 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu)) |
---|
444 | t12 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
---|
445 | t13 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2* Rg2))) + (7*b^2)/(15*q0^2*Rg2)))) |
---|
446 | t14 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
---|
447 | t15 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu))) |
---|
448 | |
---|
449 | |
---|
450 | yy = (q0^p1*(((-((b*pi)/(L*q0))) +t1/L +t2/(q0^4*Rg2^2) + 1/2*t3*t4)) + (t5*((q0^(p1 - p2)*(((-q0^(-p1))*(((b^2*pi)/(L*q0^2) +t6/L +t7/(2*C5) -t8/(C5*q0^4*Rg2^2) +t9/(q0^4*Rg2^2) -t10/(q0^5*Rg2^2) + 1/2*t11*t12)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) +t13/L +t14/(q0^4*Rg2^2) + 1/2*t15*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))) |
---|
451 | |
---|
452 | |
---|
453 | return yy |
---|
454 | end |
---|
455 | |
---|
456 | |
---|
457 | // unused functions copied from WRC Matlab code |
---|
458 | |
---|
459 | |
---|
460 | //function Vmic(phi) |
---|
461 | // Variable phi |
---|
462 | // |
---|
463 | // Variable yy |
---|
464 | // yy = 2.53*10^6*(phi)^(0.5) + 8.35*(phi)^(-2) |
---|
465 | // |
---|
466 | // return yy |
---|
467 | //end |
---|
468 | |
---|
469 | //function Scs(x) |
---|
470 | // Variable x |
---|
471 | // |
---|
472 | // Variable yy |
---|
473 | // yy = (2*bessj(1,x)/x)^2 |
---|
474 | // |
---|
475 | // return yy |
---|
476 | //end |
---|
477 | |
---|
478 | //function Srod(q,L,b) |
---|
479 | // Variable q,L,b |
---|
480 | // |
---|
481 | // Variable yy |
---|
482 | // yy = 2 * Si(q,L,b)/(q*L) - 4 * (sin(q*L/2)^2)/((q*L)^2) |
---|
483 | // |
---|
484 | // return yy |
---|
485 | //end |
---|
486 | // |
---|
487 | //function Si(q,L,b) |
---|
488 | // Variable q,L,b |
---|
489 | // |
---|
490 | // Variable yy |
---|
491 | // |
---|
492 | //// for i=1:length(q) |
---|
493 | //// y(i) = quadl('sin(x)./(x)',10e-8,q(i)*L); |
---|
494 | //// end |
---|
495 | // |
---|
496 | // return yy |
---|
497 | //end |
---|
498 | |
---|
499 | ////////////////// |
---|
500 | |
---|
501 | //wrapper to calculate the smeared model as an AAO-Struct |
---|
502 | // fills the struct and calls the ususal function with the STRUCT parameter |
---|
503 | // |
---|
504 | // used only for the dependency, not for fitting |
---|
505 | // |
---|
506 | Function fSmearedFlexExclVolCyl(coefW,yW,xW) |
---|
507 | Wave coefW,yW,xW |
---|
508 | |
---|
509 | String str = getWavesDataFolder(yW,0) |
---|
510 | String DF="root:"+str+":" |
---|
511 | |
---|
512 | WAVE resW = $(DF+str+"_res") |
---|
513 | |
---|
514 | STRUCT ResSmearAAOStruct fs |
---|
515 | WAVE fs.coefW = coefW |
---|
516 | WAVE fs.yW = yW |
---|
517 | WAVE fs.xW = xW |
---|
518 | WAVE fs.resW = resW |
---|
519 | |
---|
520 | Variable err |
---|
521 | err = SmearedFlexExclVolCyl(fs) |
---|
522 | |
---|
523 | return (0) |
---|
524 | End |
---|
525 | |
---|
526 | // this is all there is to the smeared calculation! |
---|
527 | Function SmearedFlexExclVolCyl(s) :FitFunc |
---|
528 | Struct ResSmearAAOStruct &s |
---|
529 | |
---|
530 | // the name of your unsmeared model (AAO) is the first argument |
---|
531 | Smear_Model_20(FlexExclVolCyl,s.coefW,s.xW,s.yW,s.resW) |
---|
532 | |
---|
533 | return(0) |
---|
534 | End |
---|