source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Barbell_v40.ipf @ 510

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 6.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////////
5//
6// calculates the scattering of a "barbell" shape - a cylinder with spherical end caps
7// of larger radius than the radius of the cylinder
8//
9// a double integral is used, both using Gaussian quadrature
10// routines that are now included with GaussUtils
11//
12// 76 point quadrature is necessary for both quadrature calls.
13//
14//
15// REFERENCE:
16//              H. Kaya, J. Appl. Cryst. (2004) 37, 223-230.
17//              H. Kaya and N-R deSouza, J. Appl. Cryst. (2004) 37, 508-509. (addenda and errata)
18//
19////////////////////////////////////////////////////
20
21//this macro sets up all the necessary parameters and waves that are
22//needed to calculate the model function.
23//
24Proc PlotBarbell(num,qmin,qmax)
25        Variable num=100, qmin=.001, qmax=.7
26        Prompt num "Enter number of data points for model: "
27        Prompt qmin "Enter minimum q-value (^1) for model: "
28        Prompt qmax "Enter maximum q-value (^1) for model: "
29        //
30        Make/O/D/n=(num) xwave_Barbell, ywave_Barbell
31        xwave_Barbell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
32        Make/O/D coef_Barbell = {1,20,400,40,1e-6,6.3e-6,0}                     //CH#2
33        make/o/t parameters_Barbell = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}   //CH#3
34        Edit parameters_Barbell, coef_Barbell
35       
36        Variable/G root:g_Barbell
37        g_Barbell := Barbell(coef_Barbell, ywave_Barbell, xwave_Barbell)
38        Display ywave_Barbell vs xwave_Barbell
39        ModifyGraph marker=29, msize=2, mode=4
40        ModifyGraph log=1
41        Label bottom "q (A\\S-1\\M)"
42        Label left "I(q) (cm\\S-1\\M)"
43        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
44       
45        AddModelToStrings("Barbell","coef_Barbell","parameters_Barbell","Barbell")
46//
47End
48
49
50// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
51Proc PlotSmearedBarbell(str)                                                           
52        String str
53        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
54       
55        // if any of the resolution waves are missing => abort
56        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
57                Abort
58        endif
59       
60        SetDataFolder $("root:"+str)
61       
62        // Setup parameter table for model function
63        Make/O/D smear_coef_Barbell = {1,20,400,40,1e-6,6.3e-6,0}               //CH#4
64        make/o/t smear_parameters_Barbell = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
65        Edit smear_parameters_Barbell,smear_coef_Barbell                                        //display parameters in a table
66       
67        // output smeared intensity wave, dimensions are identical to experimental QSIG values
68        // make extra copy of experimental q-values for easy plotting
69        Duplicate/O $(str+"_q") smeared_Barbell,smeared_qvals                           //
70        SetScale d,0,0,"1/cm",smeared_Barbell                                                   //
71                                       
72        Variable/G gs_Barbell=0
73        gs_Barbell := fSmearedBarbell(smear_coef_Barbell,smeared_Barbell,smeared_qvals) //this wrapper fills the STRUCT
74       
75        Display smeared_Barbell vs smeared_qvals                                                                        //
76        ModifyGraph log=1,marker=29,msize=2,mode=4
77        Label bottom "q (A\\S-1\\M)"
78        Label left "I(q) (cm\\S-1\\M)"
79        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
80       
81        SetDataFolder root:
82        AddModelToStrings("SmearedBarbell","smear_coef_Barbell","smear_parameters_Barbell","Barbell")
83End
84
85       
86
87//AAO version, uses XOP if available
88// simply calls the original single point calculation with
89// a wave assignment (this will behave nicely if given point ranges)
90Function Barbell(cw,yw,xw) : FitFunc
91        Wave cw,yw,xw
92       
93#if exists("BarbellX")
94        yw = BarbellX(cw,xw)
95#else
96        yw = fBarbell(cw,xw)
97#endif
98        return(0)
99End
100
101//
102// - a double integral - choose points wisely - 76 for both...
103//
104Function fBarbell(w,x) : FitFunc
105        Wave w
106        Variable x
107//       Input (fitting) variables are:
108        //[0] scale factor
109        //[1] cylinder radius (little r)
110        //[2] cylinder length (big L)
111        //[3] end cap radius (big R)
112        //[4] sld cylinder (A^-2)
113        //[5] sld solvent
114        //[6] incoherent background (cm^-1)
115//      give them nice names
116        Variable scale,contr,bkg,inten,sldc,slds
117        Variable len,rad,hDist,endRad
118        scale = w[0]
119        rad = w[1]
120        len = w[2]
121        endRad = w[3]
122        sldc = w[4]
123        slds = w[5]
124        bkg = w[6]
125       
126        hDist = sqrt(endRad^2-rad^2)   
127               
128        contr = sldc-slds
129       
130        Variable/G root:gDumTheta=0,root:gDumT=0
131       
132        inten = IntegrateFn76(Barbell_Outer,0,pi/2,w,x)
133       
134        inten /= pi*rad*rad*len + 2*pi*(2*endRad^3/3+endRad^2*hDist-hDist^3/3)          //divide by volume
135        inten *= 1e8            //convert to cm^-1
136        inten *= contr*contr
137        inten *= scale
138        inten += bkg
139       
140        Return (inten)
141End
142
143// outer integral
144// x is the q-value
145Function Barbell_Outer(w,x,dum)
146        Wave w
147        Variable x,dum
148       
149        Variable retVal
150        Variable scale,contr,bkg,inten,sldc,slds
151        Variable len,rad,hDist,endRad
152        scale = w[0]
153        rad = w[1]
154        len = w[2]
155        endRad = w[3]
156        sldc = w[4]
157        slds = w[5]
158        bkg = w[6]
159
160        hDist = sqrt(endRad^2-rad^2)   
161                       
162        NVAR dTheta = root:gDumTheta
163        NVAR dt = root:gDumT
164        dTheta = dum
165        retval = IntegrateFn76(Barbell_Inner,-hDist/endRad,1,w,x)
166       
167        Variable arg1,arg2
168        arg1 = x*len/2*cos(dum)
169        arg2 = x*rad*sin(dum)
170       
171        retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2
172       
173        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
174       
175        return(retVal)
176End
177
178//returns the value of the integrand of the inner integral
179Function Barbell_Inner(w,x,dum)
180        Wave w
181        Variable x,dum
182       
183        Variable retVal
184        Variable scale,contr,bkg,inten,sldc,slds
185        Variable len,rad,hDist,endRad
186        scale = w[0]
187        rad = w[1]
188        len = w[2]
189        endRad = w[3]
190        sldc = w[4]
191        slds = w[5]
192        bkg = w[6]
193       
194        NVAR dTheta = root:gDumTheta
195        NVAR dt = root:gDumT
196        dt = dum
197       
198        retVal = Barbell_integrand(w,x,dt,dTheta)
199       
200        retVal *= 4*pi*endRad^3
201       
202        return(retVal)
203End
204
205Function Barbell_integrand(w,x,tt,Theta)
206        Wave w
207        Variable x,tt,Theta
208       
209        Variable val,arg1,arg2
210        Variable scale,contr,bkg,inten,sldc,slds
211        Variable len,rad,hDist,endRad
212        scale = w[0]
213        rad = w[1]
214        len = w[2]
215        endRad = w[3]
216        sldc = w[4]
217        slds = w[5]
218        bkg = w[6]
219
220        hDist = sqrt(endRad^2-rad^2)   
221               
222        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
223        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
224       
225        val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2
226       
227        return(val)
228end
229
230//wrapper to calculate the smeared model as an AAO-Struct
231// fills the struct and calls the ususal function with the STRUCT parameter
232//
233// used only for the dependency, not for fitting
234//
235Function fSmearedBarbell(coefW,yW,xW)
236        Wave coefW,yW,xW
237       
238        String str = getWavesDataFolder(yW,0)
239        String DF="root:"+str+":"
240       
241        WAVE resW = $(DF+str+"_res")
242       
243        STRUCT ResSmearAAOStruct fs
244        WAVE fs.coefW = coefW   
245        WAVE fs.yW = yW
246        WAVE fs.xW = xW
247        WAVE fs.resW = resW
248       
249        Variable err
250        err = SmearedBarbell(fs)
251       
252        return (0)
253End
254               
255// this is all there is to the smeared calculation!
256//
257// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
258Function SmearedBarbell(s) :FitFunc
259        Struct ResSmearAAOStruct &s
260
261//      the name of your unsmeared model (AAO) is the first argument
262        Smear_Model_20(Barbell,s.coefW,s.xW,s.yW,s.resW)
263
264        return(0)
265End
Note: See TracBrowser for help on using the repository browser.