source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Cylinder_v40.ipf @ 633

Last change on this file since 633 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

File size: 7.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this function is for the form factor of a right circular cylinder with uniform scattering length density
9//
10// 06 NOV 98 SRK
11////////////////////////////////////////////////
12
13Proc PlotCylinderForm(num,qmin,qmax)
14        Variable num=128,qmin=0.001,qmax=0.7
15        Prompt num "Enter number of data points for model: "
16        Prompt qmin "Enter minimum q-value (A^-1) for model: "
17        Prompt qmax "Enter maximum q-value (A^-1) for model: "
18       
19        make/o/D/n=(num) xwave_cyl,ywave_cyl
20        xwave_cyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
21        make/o/D coef_cyl = {1.,20.,400,1e-6,6.3e-6,0.01}
22        make/o/t parameters_cyl = {"scale","radius (A)","length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
23        Edit parameters_cyl,coef_cyl
24        Variable/G root:g_cyl
25        g_cyl := CylinderForm(coef_cyl,ywave_cyl,xwave_cyl)
26//      ywave_cyl := CylinderForm(coef_cyl,xwave_cyl)
27        Display ywave_cyl vs xwave_cyl
28        ModifyGraph log=1,marker=29,msize=2,mode=4
29        Label bottom "q (A\\S-1\\M)"
30        Label left "Intensity (cm\\S-1\\M)"
31        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
32       
33        AddModelToStrings("CylinderForm","coef_cyl","parameters_cyl","cyl")
34End
35
36///////////////////////////////////////////////////////////
37// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
38Proc PlotSmearedCylinderForm(str)                                                               
39        String str
40        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
41       
42        // if any of the resolution waves are missing => abort
43        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
44                Abort
45        endif
46       
47        SetDataFolder $("root:"+str)
48       
49        // Setup parameter table for model function
50        make/o/D smear_coef_cyl = {1.,20.,400,1e-6,6.3e-6,0.01}
51        make/o/t smear_parameters_cyl = {"scale","radius (A)","length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
52        Edit smear_parameters_cyl,smear_coef_cyl
53       
54        // output smeared intensity wave, dimensions are identical to experimental QSIG values
55        // make extra copy of experimental q-values for easy plotting
56        Duplicate/O $(str+"_q") smeared_cyl,smeared_qvals
57        SetScale d,0,0,"1/cm",smeared_cyl       
58                                       
59        Variable/G gs_cyl=0
60        gs_cyl := fSmearedCylinderForm(smear_coef_cyl,smeared_cyl,smeared_qvals)        //this wrapper fills the STRUCT
61       
62        Display smeared_cyl vs smeared_qvals
63        ModifyGraph log=1,marker=29,msize=2,mode=4
64        Label bottom "q (A\\S-1\\M)"
65        Label left "Intensity (cm\\S-1\\M)"
66        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
67       
68        SetDataFolder root:
69        AddModelToStrings("SmearedCylinderForm","smear_coef_cyl","smear_parameters_cyl","cyl")
70End
71
72// AAO verison
73// even 100 points gets a 1.7x speedup from MultiThread
74// then the fit is 1.24x faster (1.25s vs 1.55s)
75Function CylinderForm(cw,yw,xw) : FitFunc
76        Wave cw,yw,xw
77
78//      Variable t1=StopMSTimer(-2)
79
80#if exists("CylinderFormX")
81//      yw = CylinderFormX(cw,xw)
82        MultiThread yw = CylinderFormX(cw,xw)
83#else
84        yw = fCylinderForm(cw,xw)
85#endif
86
87//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
88
89        return(0)
90End
91
92///////////////////////////////////////////////////////////////
93// unsmeared model calculation
94///////////////////////////
95Function fCylinderForm(w,x) : FitFunc
96        Wave w
97        Variable x
98       
99//The input variables are (and output)
100        //[0] scale
101        //[1] cylinder RADIUS (A)
102        //[2] total cylinder LENGTH (A)
103        //[3] sld cylinder (A^-2)
104        //[4] sld solvent
105        //[5] background (cm^-1)
106        Variable scale, radius,length,delrho,bkg,sldCyl,sldSolv
107        scale = w[0]
108        radius = w[1]
109        length = w[2]
110        sldCyl = w[3]
111        sldSolv = w[4]
112        bkg = w[5]
113        delrho = sldCyl-sldSolv
114
115//
116// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
117//
118
119// local variables
120        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
121        Variable answer
122        String weightStr,zStr
123       
124        weightStr = "gauss76wt"
125        zStr = "gauss76z"
126
127       
128//      if wt,z waves don't exist, create them
129// 20 Gauss points is not enough for cylinder calculation
130       
131        if (WaveExists($weightStr) == 0) // wave reference is not valid,
132                Make/D/N=76 $weightStr,$zStr
133                Wave w76 = $weightStr
134                Wave z76 = $zStr                // wave references to pass
135                Make76GaussPoints(w76,z76)     
136        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
137        else
138                if(exists(weightStr) > 1)
139                         Abort "wave name is already in use"    // execute if condition is false
140                endif
141                Wave w76 = $weightStr
142                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
143        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
144        endif
145
146
147// set up the integration
148        // end points and weights
149        nord = 76
150        va = 0
151        vb = Pi/2.0
152      halfheight = length/2.0
153
154// evaluate at Gauss points
155        // remember to index from 0,size-1
156
157        qq = x          //current x point is the q-value for evaluation
158      summ = 0.0                // initialize integral
159      ii=0
160      do
161                // Using 76 Gauss points
162                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
163                yyy = w76[ii] * cyl(qq, radius, halfheight, zi)
164                summ += yyy
165
166                ii+=1
167        while (ii<nord)                         // end of loop over quadrature points
168//   
169// calculate value of integral to return
170
171      answer = (vb-va)/2.0*summ
172     
173// Multiply by contrast^2
174        answer *= delrho*delrho
175//normalize by cylinder volume
176//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
177        vcyl=Pi*radius*radius*length
178        answer *= vcyl
179//convert to [cm-1]
180        answer *= 1.0e8
181//Scale
182        answer *= scale
183// add in the background
184        answer += bkg
185
186        Return (answer)
187       
188End             //End of function CylinderForm()
189
190///////////////////////////////////////////////////////////////
191Function cyl(qq,rr,h,theta)
192        Variable qq,rr,h,theta
193       
194// qq is the q-value for the calculation (1/A)
195// rr is the radius of the cylinder (A)
196// h is the HALF-LENGTH of the cylinder = L/2 (A)
197
198   //Local variables
199        Variable besarg,bj,retval,d1,t1,b1,t2,b2,be,si,siarg
200    besarg = qq*rr*sin(theta)
201        siarg = qq * h * cos(theta)
202       
203    bj =bessJ(1,besarg)
204
205//* Computing 2nd power */
206    d1 = sin(siarg)
207    t1 = d1 * d1
208//* Computing 2nd power */
209    d1 = bj
210    t2 = d1 * d1 * 4.0 * sin(theta)
211//* Computing 2nd power */
212    d1 = siarg
213    b1 = d1 * d1
214//* Computing 2nd power */
215    d1 = qq * rr * sin(theta)
216    b2 = d1 * d1
217   
218    if(besarg == 0)
219        be = sin(theta)
220    else
221        be = t2/b2
222    endif
223   
224    if(siarg == 0)
225        si = 1
226    else
227        si = t1/b1
228    endif
229   
230    retval = be*si
231
232    return retval
233   
234End     //Function cyl()
235
236// this is all there is to the smeared calculation!
237Function SmearedCylinderForm(s) :FitFunc
238        Struct ResSmearAAOStruct &s
239
240////the name of your unsmeared model is the first argument
241        Smear_Model_20(CylinderForm,s.coefW,s.xW,s.yW,s.resW)
242
243        return(0)
244End
245
246
247//wrapper to calculate the smeared model as an AAO-Struct
248// fills the struct and calls the ususal function with the STRUCT parameter
249//
250// used only for the dependency, not for fitting
251//
252Function fSmearedCylinderForm(coefW,yW,xW)
253        Wave coefW,yW,xW
254       
255        String str = getWavesDataFolder(yW,0)
256        String DF="root:"+str+":"
257       
258        WAVE resW = $(DF+str+"_res")
259       
260        STRUCT ResSmearAAOStruct fs
261        WAVE fs.coefW = coefW   
262        WAVE fs.yW = yW
263        WAVE fs.xW = xW
264        WAVE fs.resW = resW
265       
266        Variable err
267        err = SmearedCylinderForm(fs)
268       
269        return (0)
270End
Note: See TracBrowser for help on using the repository browser.