source: sans/Dev/trunk/NCNR_User_Procedures/SANS_Analysis/Models/Cylinder_v40.ipf @ 325

Last change on this file since 325 was 325, checked in by ajj, 15 years ago

Adding SANS Analysis Models to new Dev tree

File size: 6.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
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 (^-1) for model: "
17        Prompt qmax "Enter maximum q-value (^-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 (\\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","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 (\\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","cyl")
70End
71
72// AAO verison
73Function CylinderForm(cw,yw,xw) : FitFunc
74        Wave cw,yw,xw
75
76#if exists("CylinderFormX")
77        yw = CylinderFormX(cw,xw)
78#else
79        yw = fCylinderForm(cw,xw)
80#endif
81        return(0)
82End
83
84///////////////////////////////////////////////////////////////
85// unsmeared model calculation
86///////////////////////////
87Function fCylinderForm(w,x) : FitFunc
88        Wave w
89        Variable x
90       
91//The input variables are (and output)
92        //[0] scale
93        //[1] cylinder RADIUS (A)
94        //[2] total cylinder LENGTH (A)
95        //[3] sld cylinder (A^-2)
96        //[4] sld solvent
97        //[5] background (cm^-1)
98        Variable scale, radius,length,delrho,bkg,sldCyl,sldSolv
99        scale = w[0]
100        radius = w[1]
101        length = w[2]
102        sldCyl = w[3]
103        sldSolv = w[4]
104        bkg = w[5]
105        delrho = sldCyl-sldSolv
106
107//
108// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
109//
110
111// local variables
112        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
113        Variable answer
114        String weightStr,zStr
115       
116        weightStr = "gauss76wt"
117        zStr = "gauss76z"
118
119       
120//      if wt,z waves don't exist, create them
121// 20 Gauss points is not enough for cylinder calculation
122       
123        if (WaveExists($weightStr) == 0) // wave reference is not valid,
124                Make/D/N=76 $weightStr,$zStr
125                Wave w76 = $weightStr
126                Wave z76 = $zStr                // wave references to pass
127                Make76GaussPoints(w76,z76)     
128        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
129        else
130                if(exists(weightStr) > 1)
131                         Abort "wave name is already in use"    // execute if condition is false
132                endif
133                Wave w76 = $weightStr
134                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
135        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
136        endif
137
138
139// set up the integration
140        // end points and weights
141        nord = 76
142        va = 0
143        vb = Pi/2.0
144      halfheight = length/2.0
145
146// evaluate at Gauss points
147        // remember to index from 0,size-1
148
149        qq = x          //current x point is the q-value for evaluation
150      summ = 0.0                // initialize integral
151      ii=0
152      do
153                // Using 76 Gauss points
154                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
155                yyy = w76[ii] * cyl(qq, radius, halfheight, zi)
156                summ += yyy
157
158                ii+=1
159        while (ii<nord)                         // end of loop over quadrature points
160//   
161// calculate value of integral to return
162
163      answer = (vb-va)/2.0*summ
164     
165// Multiply by contrast^2
166        answer *= delrho*delrho
167//normalize by cylinder volume
168//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
169        vcyl=Pi*radius*radius*length
170        answer *= vcyl
171//convert to [cm-1]
172        answer *= 1.0e8
173//Scale
174        answer *= scale
175// add in the background
176        answer += bkg
177
178        Return (answer)
179       
180End             //End of function CylinderForm()
181
182///////////////////////////////////////////////////////////////
183Function cyl(qq,rr,h,theta)
184        Variable qq,rr,h,theta
185       
186// qq is the q-value for the calculation (1/A)
187// rr is the radius of the cylinder (A)
188// h is the HALF-LENGTH of the cylinder = L/2 (A)
189
190   //Local variables
191        Variable besarg,bj,retval,d1,t1,b1,t2,b2
192    besarg = qq*rr*sin(theta)
193
194    bj =bessJ(1,besarg)
195
196//* Computing 2nd power */
197    d1 = sin(qq * h * cos(theta))
198    t1 = d1 * d1
199//* Computing 2nd power */
200    d1 = bj
201    t2 = d1 * d1 * 4.0 * sin(theta)
202//* Computing 2nd power */
203    d1 = qq * h * cos(theta)
204    b1 = d1 * d1
205//* Computing 2nd power */
206    d1 = qq * rr * sin(theta)
207    b2 = d1 * d1
208    retval = t1 * t2 / b1 / b2
209
210    return retval
211   
212End     //Function cyl()
213
214// this is all there is to the smeared calculation!
215Function SmearedCylinderForm(s) :FitFunc
216        Struct ResSmearAAOStruct &s
217
218////the name of your unsmeared model is the first argument
219        Smear_Model_20(CylinderForm,s.coefW,s.xW,s.yW,s.resW)
220
221        return(0)
222End
223
224
225//wrapper to calculate the smeared model as an AAO-Struct
226// fills the struct and calls the ususal function with the STRUCT parameter
227//
228// used only for the dependency, not for fitting
229//
230Function fSmearedCylinderForm(coefW,yW,xW)
231        Wave coefW,yW,xW
232       
233        String str = getWavesDataFolder(yW,0)
234        String DF="root:"+str+":"
235       
236        WAVE resW = $(DF+str+"_res")
237       
238        STRUCT ResSmearAAOStruct fs
239        WAVE fs.coefW = coefW   
240        WAVE fs.yW = yW
241        WAVE fs.xW = xW
242        WAVE fs.resW = resW
243       
244        Variable err
245        err = SmearedCylinderForm(fs)
246       
247        return (0)
248End
Note: See TracBrowser for help on using the repository browser.