source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/CylinderForm.ipf @ 151

Last change on this file since 151 was 151, checked in by srkline, 15 years ago

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

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