source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Spherocylinder_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: 7.0 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 spherocylinder, that is a cylinder with spherical end caps
7// where the radius of the end caps is the same as 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 PlotSpherocylinder(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_SphCyl, ywave_SphCyl
31        xwave_SphCyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
32        Make/O/D coef_SphCyl = {1,20,400,1e-6,6.3e-6,0}                 //CH#2
33        make/o/t parameters_SphCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"} //CH#3
34        Edit parameters_SphCyl, coef_SphCyl
35       
36        Variable/G root:g_SphCyl
37        g_SphCyl := Spherocylinder(coef_SphCyl, ywave_SphCyl, xwave_SphCyl)
38        Display ywave_SphCyl vs xwave_SphCyl
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("Spherocylinder","coef_SphCyl","parameters_SphCyl","SphCyl")
46//
47End
48
49
50// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
51Proc PlotSmearedSpherocylinder(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_SphCyl = {1,20,400,1e-6,6.3e-6,0}           //CH#4
64        make/o/t smear_parameters_SphCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
65        Edit smear_parameters_SphCyl,smear_coef_SphCyl                                  //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_SphCyl,smeared_qvals                            //
70        SetScale d,0,0,"1/cm",smeared_SphCyl                                                    //
71                                       
72        Variable/G gs_SphCyl=0
73        gs_SphCyl := fSmearedSpherocylinder(smear_coef_SphCyl,smeared_SphCyl,smeared_qvals)     //this wrapper fills the STRUCT
74       
75        Display smeared_SphCyl 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("SmearedSpherocylinder","smear_coef_SphCyl","smear_parameters_SphCyl","SphCyl")
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 Spherocylinder(cw,yw,xw) : FitFunc
91        Wave cw,yw,xw
92       
93#if exists("SpherocylinderX")
94        yw = SpherocylinderX(cw,xw)
95#else
96        yw = fSpherocylinder(cw,xw)
97#endif
98        return(0)
99End
100
101//
102// - a double integral - choose points wisely - 76 for both...
103//
104Function fSpherocylinder(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[3]
123        slds = w[4]
124        bkg = w[5]
125       
126        Make/O/D/N=7 SphCyl_tmp
127        SphCyl_tmp[0] = w[0]
128        SphCyl_tmp[1] = w[1]
129        SphCyl_tmp[2] = w[2]
130        SphCyl_tmp[3] = w[1]            //end radius is same as cylinder radius
131        SphCyl_tmp[4] = w[3]
132        SphCyl_tmp[5] = w[4]
133        SphCyl_tmp[6] = w[5]
134       
135        hDist = 0               //by definition
136               
137        contr = sldc-slds
138       
139        Variable/G root:gDumTheta=0,root:gDumT=0
140       
141        inten = IntegrateFn76(SphCyl_Outer,0,pi/2,SphCyl_tmp,x)
142       
143        inten /= pi*rad*rad*len + pi*4*endRad^3/3               //divide by volume
144        inten *= 1e8            //convert to cm^-1
145        inten *= contr*contr
146        inten *= scale
147        inten += bkg
148       
149        Return (inten)
150End
151
152// outer integral
153// x is the q-value
154Function SphCyl_Outer(w,x,dum)
155        Wave w
156        Variable x,dum
157       
158        Variable retVal
159        Variable scale,contr,bkg,inten,sldc,slds
160        Variable len,rad,hDist,endRad
161        scale = w[0]
162        rad = w[1]
163        len = w[2]
164        endRad = w[3]
165        sldc = w[4]
166        slds = w[5]
167        bkg = w[6]
168
169        hDist = 0
170                               
171        NVAR dTheta = root:gDumTheta
172        NVAR dt = root:gDumT
173        dTheta = dum
174        retval = IntegrateFn76(SphCyl_Inner,-hDist/endRad,1,w,x)
175       
176        Variable arg1,arg2
177        arg1 = x*len/2*cos(dum)
178        arg2 = x*rad*sin(dum)
179       
180        retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2
181       
182        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
183       
184        return(retVal)
185End
186
187//returns the value of the integrand of the inner integral
188Function SphCyl_Inner(w,x,dum)
189        Wave w
190        Variable x,dum
191       
192        Variable retVal
193        Variable scale,contr,bkg,inten,sldc,slds
194        Variable len,rad,hDist,endRad
195        scale = w[0]
196        rad = w[1]
197        len = w[2]
198        endRad = w[3]
199        sldc = w[4]
200        slds = w[5]
201        bkg = w[6]
202       
203        NVAR dTheta = root:gDumTheta
204        NVAR dt = root:gDumT
205        dt = dum
206       
207        retVal = SphCyl(w,x,dt,dTheta)
208       
209        retVal *= 4*pi*endRad^3
210       
211        return(retVal)
212End
213
214Function SphCyl(w,x,tt,Theta)
215        Wave w
216        Variable x,tt,Theta
217       
218        Variable val,arg1,arg2
219        Variable scale,contr,bkg,inten,sldc,slds
220        Variable len,rad,hDist,endRad
221        scale = w[0]
222        rad = w[1]
223        len = w[2]
224        endRad = w[3]
225        sldc = w[4]
226        slds = w[5]
227        bkg = w[6]
228       
229        hDist = 0
230               
231        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
232        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
233       
234        val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2
235       
236        return(val)
237end
238
239//wrapper to calculate the smeared model as an AAO-Struct
240// fills the struct and calls the ususal function with the STRUCT parameter
241//
242// used only for the dependency, not for fitting
243//
244Function fSmearedSpherocylinder(coefW,yW,xW)
245        Wave coefW,yW,xW
246       
247        String str = getWavesDataFolder(yW,0)
248        String DF="root:"+str+":"
249       
250        WAVE resW = $(DF+str+"_res")
251       
252        STRUCT ResSmearAAOStruct fs
253        WAVE fs.coefW = coefW   
254        WAVE fs.yW = yW
255        WAVE fs.xW = xW
256        WAVE fs.resW = resW
257       
258        Variable err
259        err = SmearedSpherocylinder(fs)
260       
261        return (0)
262End
263               
264// this is all there is to the smeared calculation!
265//
266// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
267Function SmearedSpherocylinder(s) :FitFunc
268        Struct ResSmearAAOStruct &s
269
270//      the name of your unsmeared model (AAO) is the first argument
271        Smear_Model_20(Spherocylinder,s.coefW,s.xW,s.yW,s.resW)
272
273        return(0)
274End
Note: See TracBrowser for help on using the repository browser.