source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/CylinderForm.ipf @ 381

Last change on this file since 381 was 381, checked in by srkline, 14 years ago

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 6.0 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 (A^-1) for model: "
16        Prompt qmax "Enter maximum q-value (A^-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 (A\\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 (A\\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// unsmeared model calculation
61///////////////////////////
62Function CylinderForm(w,x) : FitFunc
63        Wave w
64        Variable x
65
66//#if exists("CylinderFormX")
67//      return CylinderFormX(w,x)
68//#endif
69
70//      if(exists("CylinderFormX")==3)
71//              FUNCREF SANSModel_proto func=$"CylinderFormX"
72//              return func(w,x)                        //defined this way to hide the functionX name from the compiler
73//      endif
74       
75//The input variables are (and output)
76        //[0] scale
77        //[1] cylinder RADIUS (A)
78        //[2] total cylinder LENGTH (A)
79        //[3] contrast (A^-2)
80        //[4] background (cm^-1)
81        Variable scale, radius,length,delrho,bkg
82        scale = w[0]
83        radius = w[1]
84        length = w[2]
85        delrho = w[3]
86        bkg = w[4]
87//
88// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
89//
90
91// local variables
92        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
93        Variable answer
94        String weightStr,zStr
95       
96        weightStr = "gauss76wt"
97        zStr = "gauss76z"
98
99       
100//      if wt,z waves don't exist, create them
101// 20 Gauss points is not enough for cylinder calculation
102       
103        if (WaveExists($weightStr) == 0) // wave reference is not valid,
104                Make/D/N=76 $weightStr,$zStr
105                Wave w76 = $weightStr
106                Wave z76 = $zStr                // wave references to pass
107                Make76GaussPoints(w76,z76)     
108        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
109        else
110                if(exists(weightStr) > 1)
111                         Abort "wave name is already in use"    // execute if condition is false
112                endif
113                Wave w76 = $weightStr
114                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
115        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
116        endif
117
118
119// set up the integration
120        // end points and weights
121        nord = 76
122        va = 0
123        vb = Pi/2.0
124      halfheight = length/2.0
125
126// evaluate at Gauss points
127        // remember to index from 0,size-1
128
129        qq = x          //current x point is the q-value for evaluation
130      summ = 0.0                // initialize integral
131      ii=0
132      do
133                // Using 76 Gauss points
134                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
135                yyy = w76[ii] * cyl(qq, radius, halfheight, zi)
136                summ += yyy
137
138                ii+=1
139        while (ii<nord)                         // end of loop over quadrature points
140//   
141// calculate value of integral to return
142
143      answer = (vb-va)/2.0*summ
144     
145// Multiply by contrast^2
146        answer *= delrho*delrho
147//normalize by cylinder volume
148//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
149        vcyl=Pi*radius*radius*length
150        answer *= vcyl
151//convert to [cm-1]
152        answer *= 1.0e8
153//Scale
154        answer *= scale
155// add in the background
156        answer += bkg
157
158        Return (answer)
159       
160End             //End of function CylinderForm()
161
162///////////////////////////////////////////////////////////////
163Function cyl(qq,rr,h,theta)
164        Variable qq,rr,h,theta
165       
166// qq is the q-value for the calculation (1/A)
167// rr is the radius of the cylinder (A)
168// h is the HALF-LENGTH of the cylinder = L/2 (A)
169
170   //Local variables
171        Variable besarg,bj,retval,d1,t1,b1,t2,b2
172    besarg = qq*rr*sin(theta)
173
174    bj =bessJ(1,besarg)
175
176//* Computing 2nd power */
177    d1 = sin(qq * h * cos(theta))
178    t1 = d1 * d1
179//* Computing 2nd power */
180    d1 = bj
181    t2 = d1 * d1 * 4.0 * sin(theta)
182//* Computing 2nd power */
183    d1 = qq * h * cos(theta)
184    b1 = d1 * d1
185//* Computing 2nd power */
186    d1 = qq * rr * sin(theta)
187    b2 = d1 * d1
188    retval = t1 * t2 / b1 / b2
189
190    return retval
191   
192End     //Function cyl()
193
194// this is all there is to the smeared calculation!
195Function SmearedCylinderForm(w,x) :FitFunc
196        Wave w
197        Variable x
198       
199        Variable ans
200        SVAR sq = gSig_Q
201        SVAR qb = gQ_bar
202        SVAR sh = gShadow
203        SVAR gQ = gQVals
204       
205        //the name of your unsmeared model is the first argument
206//      if(exists("CylinderFormX") == 3)
207//              ans = Smear_Model_20($"CylinderFormX",$sq,$qb,$sh,$gQ,w,x)
208//      else
209                ans = Smear_Model_20(CylinderForm,$sq,$qb,$sh,$gQ,w,x)
210//      endif
211
212        return(ans)
213End
214
Note: See TracBrowser for help on using the repository browser.