source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Cylinder_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.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 (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
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.