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

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

initial checkin of Analysis v.3.00 files

File size: 5.4 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 an ellipsoid of rotation with uniform scattering length density
8//
9// 06 NOV 98 SRK
10////////////////////////////////////////////////
11
12Proc PlotEllipsoidForm(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_eor,ywave_eor
19        xwave_eor =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
20        Make/O/D coef_eor = {1.,20.,400,3.0e-6,0.01}
21        make/o/t parameters_eor = {"scale","R a (rotation axis) (A)","R b (A)","Contrast (A^-2)","incoh. bkg (cm^-1)"}
22        Edit parameters_eor,coef_eor
23        ywave_eor := EllipsoidForm(coef_eor,xwave_eor)
24        Display ywave_eor vs xwave_eor
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 PlotSmearedEllipsoidForm()
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_eor = {1.,20.,400,3.0e-6,0.01}
43        make/o/t smear_parameters_eor = {"scale","R a (rotation axis) (A)","R b (A)","Contrast (A^-2)","incoh. bkg (cm^-1)"}
44        Edit smear_parameters_eor,smear_coef_eor
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_eor,smeared_qvals
49        SetScale d,0,0,"1/cm",smeared_eor       
50
51        smeared_eor := SmearedEllipsoidForm(smear_coef_eor,$gQvals)
52        Display smeared_eor 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 EllipsoidForm(w,x) : FitFunc
63        Wave w
64        Variable x
65
66//The input variables are (and output)
67        //[0] scale
68        //[1] Axis of rotation
69        //[2] two equal radii
70        //[3] contrast (A^-2)
71        //[4] background (cm^-1)
72        Variable scale, ra,vra,delrho,bkg
73        scale = w[0]
74        vra = w[1]
75        ra = w[2]
76        delrho = w[3]
77        bkg = w[4]
78        //if vra < ra, OBLATE
79        //if vra > ra, PROLATE
80//
81// the OUTPUT form factor is <f^2>/Vell [cm-1]
82//
83
84// local variables
85        Variable nord,ii,va,vb,contr,vell,nden,summ,yyy,zi,qq,halfheight
86        Variable answer
87        String weightStr,zStr
88       
89        weightStr = "gauss76wt"
90        zStr = "gauss76z"
91
92       
93//      if wt,z waves don't exist, create them
94// 20 Gauss points is not enough for Ellipsoid calculation
95       
96        if (WaveExists($weightStr) == 0) // wave reference is not valid,
97                Make/D/N=76 $weightStr,$zStr
98                Wave w76 = $weightStr
99                Wave z76 = $zStr                // wave references to pass
100                Make76GaussPoints(w76,z76)     
101        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
102        else
103                if(exists(weightStr) > 1)
104                         Abort "wave name is already in use"    // execute if condition is false
105                endif
106                Wave w76 = $weightStr
107                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
108        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
109        endif
110
111
112// set up the integration
113        // end points and weights
114        nord = 76
115        va = 0
116        vb = 1
117
118// evaluate at Gauss points
119        // remember to index from 0,size-1
120
121        qq = x          //current x point is the q-value for evaluation
122      summ = 0.0                // initialize integral
123      ii=0
124      do
125                // Using 76 Gauss points
126                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
127                yyy = w76[ii] * eor(qq, ra,vra, zi)
128                summ += yyy
129
130                ii+=1
131        while (ii<nord)                         // end of loop over quadrature points
132//   
133// calculate value of integral to return
134
135      answer = (vb-va)/2.0*summ
136     
137// Multiply by contrast^2
138        answer *= delrho*delrho
139//normalize by Ellipsoid volume
140//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vell
141        vell=4*Pi/3*ra*ra*vra
142        answer *= vell
143//convert to [cm-1]
144        answer *= 1.0e8
145//Scale
146        answer *= scale
147// add in the background
148        answer += bkg
149
150        Return (answer)
151       
152End             //End of function EllipsoidForm()
153
154///////////////////////////////////////////////////////////////
155Function eor(qq,ra,vra,theta)
156        Variable qq,ra,vra,theta
157       
158// qq is the q-value for the calculation (1/A)
159// ra are the like radius of the Ellipsoid (A)
160// vra is the unlike semiaxis of the Ellipsoid = L/2 (A) (= the rotation axis!)
161// theta is the dummy variable of the integration
162
163   //Local variables
164        Variable retval,arg,t1, nu
165       
166        nu = vra/ra
167        arg = qq*ra*sqrt(1+theta^2*(nu^2-1))
168       
169        retval = 9*((sin(arg)-arg*cos(arg))/arg^3)^2
170       
171    return retval
172   
173End     //Function eor()
174
175// this is all there is to the smeared calculation!
176Function SmearedEllipsoidForm(w,x) :FitFunc
177        Wave w
178        Variable x
179       
180        Variable ans
181        SVAR sq = gSig_Q
182        SVAR qb = gQ_bar
183        SVAR sh = gShadow
184        SVAR gQ = gQVals
185       
186        //the name of your unsmeared model is the first argument
187        ans = Smear_Model_20(EllipsoidForm,$sq,$qb,$sh,$gQ,w,x)
188
189        return(ans)
190End
Note: See TracBrowser for help on using the repository browser.