source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/CylinderForm.ipf @ 42

Last change on this file since 42 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

File size: 5.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
5// Adopting these into the experiment will insure that they are always present
6////////////////////////////////////////////////
7// this function is for the form factor of a right circular cylinder with uniform scattering length density
8//
9// 06 NOV 98 SRK
10////////////////////////////////////////////////
11
12Proc PlotCylinderForm(num,qmin,qmax)
13        Variable num=128,qmin=0.001,qmax=0.7
14        Prompt num "Enter number of data points for model: "
15        Prompt qmin "Enter minimum q-value (^-1) for model: "
16        Prompt qmax "Enter maximum q-value (^-1) for model: "
17       
18        make/o/D/n=(num) xwave_cyl,ywave_cyl
19        xwave_cyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
20        make/o/D coef_cyl = {1.,20.,400,3.0e-6,0.01}
21        make/o/t parameters_cyl = {"scale","radius (A)","length (A)","contrast (A^-2)","incoh. bkg (cm^-1)"}
22        Edit parameters_cyl,coef_cyl
23        ywave_cyl := CylinderForm(coef_cyl,xwave_cyl)
24        Display ywave_cyl vs xwave_cyl
25        ModifyGraph log=1,marker=29,msize=2,mode=4
26        Label bottom "q (\\S-1\\M)"
27        Label left "Intensity (cm\\S-1\\M)"
28        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
29End
30///////////////////////////////////////////////////////////
31
32Proc PlotSmearedCylinderForm() 
33        //no input parameters necessary, it MUST use the experimental q-values
34        // from the experimental data read in from an AVE/QSIG data file
35       
36        // if no gQvals wave, data must not have been loaded => abort
37        if(ResolutionWavesMissing())
38                Abort
39        endif
40       
41        // Setup parameter table for model function
42        make/o/D smear_coef_cyl = {1.,20.,400,3.0e-6,0.01}
43        make/o/t smear_parameters_cyl = {"scale","radius (A)","length (A)","contrast (A^-2)","incoh. bkg (cm^-1)"}
44        Edit smear_parameters_cyl,smear_coef_cyl
45       
46        // output smeared intensity wave, dimensions are identical to experimental QSIG values
47        // make extra copy of experimental q-values for easy plotting
48        Duplicate/O $gQvals smeared_cyl,smeared_qvals
49        SetScale d,0,0,"1/cm",smeared_cyl       
50
51        smeared_cyl := SmearedCylinderForm(smear_coef_cyl,$gQvals)
52        Display smeared_cyl vs smeared_qvals
53        ModifyGraph log=1,marker=29,msize=2,mode=4
54        Label bottom "q (\\S-1\\M)"
55        Label left "Intensity (cm\\S-1\\M)"
56        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
57End
58
59
60///////////////////////////////////////////////////////////////
61// unsmeared model calculation
62///////////////////////////
63Function CylinderForm(w,x) : FitFunc
64        Wave w
65        Variable x
66
67        String funcStr = SelectString(exists("CylinderFormX")==3,"", "CylinderFormX")
68        if(strlen(funcStr) > 0)
69                FUNCREF SANSModel_proto func=$funcStr
70                return func(w,x)
71        endif
72       
73//The input variables are (and output)
74        //[0] scale
75        //[1] cylinder RADIUS (A)
76        //[2] total cylinder LENGTH (A)
77        //[3] contrast (A^-2)
78        //[4] background (cm^-1)
79        Variable scale, radius,length,delrho,bkg
80        scale = w[0]
81        radius = w[1]
82        length = w[2]
83        delrho = w[3]
84        bkg = w[4]
85//
86// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
87//
88
89// local variables
90        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
91        Variable answer
92        String weightStr,zStr
93       
94        weightStr = "gauss76wt"
95        zStr = "gauss76z"
96
97       
98//      if wt,z waves don't exist, create them
99// 20 Gauss points is not enough for cylinder calculation
100       
101        if (WaveExists($weightStr) == 0) // wave reference is not valid,
102                Make/D/N=76 $weightStr,$zStr
103                Wave w76 = $weightStr
104                Wave z76 = $zStr                // wave references to pass
105                Make76GaussPoints(w76,z76)     
106        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
107        else
108                if(exists(weightStr) > 1)
109                         Abort "wave name is already in use"    // execute if condition is false
110                endif
111                Wave w76 = $weightStr
112                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
113        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
114        endif
115
116
117// set up the integration
118        // end points and weights
119        nord = 76
120        va = 0
121        vb = Pi/2.0
122      halfheight = length/2.0
123
124// evaluate at Gauss points
125        // remember to index from 0,size-1
126
127        qq = x          //current x point is the q-value for evaluation
128      summ = 0.0                // initialize integral
129      ii=0
130      do
131                // Using 76 Gauss points
132                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
133                yyy = w76[ii] * cyl(qq, radius, halfheight, zi)
134                summ += yyy
135
136                ii+=1
137        while (ii<nord)                         // end of loop over quadrature points
138//   
139// calculate value of integral to return
140
141      answer = (vb-va)/2.0*summ
142     
143// Multiply by contrast^2
144        answer *= delrho*delrho
145//normalize by cylinder volume
146//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
147        vcyl=Pi*radius*radius*length
148        answer *= vcyl
149//convert to [cm-1]
150        answer *= 1.0e8
151//Scale
152        answer *= scale
153// add in the background
154        answer += bkg
155
156        Return (answer)
157       
158End             //End of function CylinderForm()
159
160///////////////////////////////////////////////////////////////
161Function cyl(qq,rr,h,theta)
162        Variable qq,rr,h,theta
163       
164// qq is the q-value for the calculation (1/A)
165// rr is the radius of the cylinder (A)
166// h is the HALF-LENGTH of the cylinder = L/2 (A)
167
168   //Local variables
169        Variable besarg,bj,retval,d1,t1,b1,t2,b2
170    besarg = qq*rr*sin(theta)
171
172    bj =bessJ(1,besarg)
173
174//* Computing 2nd power */
175    d1 = sin(qq * h * cos(theta))
176    t1 = d1 * d1
177//* Computing 2nd power */
178    d1 = bj
179    t2 = d1 * d1 * 4.0 * sin(theta)
180//* Computing 2nd power */
181    d1 = qq * h * cos(theta)
182    b1 = d1 * d1
183//* Computing 2nd power */
184    d1 = qq * rr * sin(theta)
185    b2 = d1 * d1
186    retval = t1 * t2 / b1 / b2
187
188    return retval
189   
190End     //Function cyl()
191
192// this is all there is to the smeared calculation!
193Function SmearedCylinderForm(w,x) :FitFunc
194        Wave w
195        Variable x
196       
197        Variable ans
198        SVAR sq = gSig_Q
199        SVAR qb = gQ_bar
200        SVAR sh = gShadow
201        SVAR gQ = gQVals
202       
203        //the name of your unsmeared model is the first argument
204        if(exists("CylinderFormX") == 3)
205                ans = Smear_Model_20($"CylinderFormX",$sq,$qb,$sh,$gQ,w,x)
206        else
207                ans = Smear_Model_20(CylinderForm,$sq,$qb,$sh,$gQ,w,x)
208        endif
209
210        return(ans)
211End
212
Note: See TracBrowser for help on using the repository browser.