source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexibleCylinder_v40.ipf @ 570

Revision 570, 15.1 KB checked in by srkline, 5 years ago (diff)

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
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//
21Proc 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 (A^-1) for model: "
25        Prompt qmax "Enter maximum q-value (A^-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 (A\\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","parameters_fle","fle")
42End
43///////////////////////////////////////////////////////////
44
45// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
46Proc 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 (A\\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","smear_parameters_fle","fle")
78End
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)
84Function 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)
93End
94
95//
96Function 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)
136End
137
138//////////////////WRC corrected code below
139// main function
140function 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)
186end
187
188//WR named this w (too generic)
189Function 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)
198end
199
200//
201function 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
208end
209
210// was named u
211function 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
217end
218
219
220
221//
222function 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
229end
230
231//
232function 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
239end
240
241//
242function 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
249end
250
251//
252function 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
259end
260
261//
262function 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
269end
270
271//
272function 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
279end
280
281//
282function 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
289end
290
291//
292function 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
306end
307
308// this must be WR modified version
309function 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
332end
333
334
335
336// these are the messy ones
337function 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
350end
351
352//
353function 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
365end
366
367// this one will be lots of trouble
368function 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
404end
405
406//need to define this on my own
407Function sech_WR(x)
408        variable x
409       
410        return(1/cosh(x))
411end
412//
413function 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
454end
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//
506Function 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)
524End
525
526// this is all there is to the smeared calculation!
527Function 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)
534End
Note: See TracBrowser for help on using the repository browser.