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

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

modified to switch between the XOP (if it exists) or the
Igor code if the XOP is not present

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