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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

File size: 6.3 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 an ellipsoid of rotation with uniform scattering length density
9//
10// 06 NOV 98 SRK
11////////////////////////////////////////////////
12
13Proc PlotEllipsoidForm(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_eor,ywave_eor
20        xwave_eor =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
21        Make/O/D coef_eor = {1.,20.,400,3.0e-6,0.01}
22        make/o/t parameters_eor = {"scale","R a (rotation axis) (A)","R b (A)","Contrast (A^-2)","incoh. bkg (cm^-1)"}
23        Edit parameters_eor,coef_eor
24        Variable/G root:g_eor
25        g_eor := EllipsoidForm(coef_eor,ywave_eor,xwave_eor)
26//      ywave_eor := EllipsoidForm(coef_eor,xwave_eor)
27        Display ywave_eor vs xwave_eor
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 PlotSmearedEllipsoidForm(str)     
37        String str
38        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
39       
40        // if no gQvals wave, data must not have been loaded => abort
41        if(ResolutionWavesMissingDF(str))
42                Abort
43        endif
44       
45        SetDataFolder $("root:"+str)
46       
47        // Setup parameter table for model function
48        Make/O/D smear_coef_eor = {1.,20.,400,3.0e-6,0.01}
49        make/o/t smear_parameters_eor = {"scale","R a (rotation axis) (A)","R b (A)","Contrast (A^-2)","incoh. bkg (cm^-1)"}
50        Edit smear_parameters_eor,smear_coef_eor
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_eor,smeared_qvals
55        SetScale d,0,0,"1/cm",smeared_eor       
56
57        Variable/G gs_eor=0
58        gs_eor := fSmearedEllipsoidForm(smear_coef_eor,smeared_eor,smeared_qvals)
59       
60        Display smeared_eor 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 version
70Function EllipsoidForm(cw,yw,xw) : FitFunc
71        Wave cw,yw,xw
72       
73#if exists("EllipsoidFormX")
74        yw = EllipsoidFormX(cw,xw)
75#else
76        yw = fEllipsoidForm(cw,xw)
77#endif
78        return(0)
79End
80
81///////////////////////////////////////////////////////////////
82// unsmeared model calculation
83///////////////////////////
84Function fEllipsoidForm(w,x) : FitFunc
85        Wave w
86        Variable x
87
88//The input variables are (and output)
89        //[0] scale
90        //[1] Axis of rotation
91        //[2] two equal radii
92        //[3] contrast (A^-2)
93        //[4] background (cm^-1)
94        Variable scale, ra,vra,delrho,bkg
95        scale = w[0]
96        vra = w[1]
97        ra = w[2]
98        delrho = w[3]
99        bkg = w[4]
100        //if vra < ra, OBLATE
101        //if vra > ra, PROLATE
102//
103// the OUTPUT form factor is <f^2>/Vell [cm-1]
104//
105
106// local variables
107        Variable nord,ii,va,vb,contr,vell,nden,summ,yyy,zi,qq,halfheight
108        Variable answer
109        String weightStr,zStr
110       
111        weightStr = "gauss76wt"
112        zStr = "gauss76z"
113
114       
115//      if wt,z waves don't exist, create them
116// 20 Gauss points is not enough for Ellipsoid calculation
117       
118        if (WaveExists($weightStr) == 0) // wave reference is not valid,
119                Make/D/N=76 $weightStr,$zStr
120                Wave w76 = $weightStr
121                Wave z76 = $zStr                // wave references to pass
122                Make76GaussPoints(w76,z76)     
123        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
124        else
125                if(exists(weightStr) > 1)
126                         Abort "wave name is already in use"    // execute if condition is false
127                endif
128                Wave w76 = $weightStr
129                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
130        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
131        endif
132
133
134// set up the integration
135        // end points and weights
136        nord = 76
137        va = 0
138        vb = 1
139
140// evaluate at Gauss points
141        // remember to index from 0,size-1
142
143        qq = x          //current x point is the q-value for evaluation
144      summ = 0.0                // initialize integral
145      ii=0
146      do
147                // Using 76 Gauss points
148                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
149                yyy = w76[ii] * eor(qq, ra,vra, zi)
150                summ += yyy
151
152                ii+=1
153        while (ii<nord)                         // end of loop over quadrature points
154//   
155// calculate value of integral to return
156
157      answer = (vb-va)/2.0*summ
158     
159// Multiply by contrast^2
160        answer *= delrho*delrho
161//normalize by Ellipsoid volume
162//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vell
163        vell=4*Pi/3*ra*ra*vra
164        answer *= vell
165//convert to [cm-1]
166        answer *= 1.0e8
167//Scale
168        answer *= scale
169// add in the background
170        answer += bkg
171
172        Return (answer)
173       
174End             //End of function EllipsoidForm()
175
176///////////////////////////////////////////////////////////////
177Function eor(qq,ra,vra,theta)
178        Variable qq,ra,vra,theta
179       
180// qq is the q-value for the calculation (1/A)
181// ra are the like radius of the Ellipsoid (A)
182// vra is the unlike semiaxis of the Ellipsoid = L/2 (A) (= the rotation axis!)
183// theta is the dummy variable of the integration
184
185   //Local variables
186        Variable retval,arg,t1, nu
187       
188        nu = vra/ra
189        arg = qq*ra*sqrt(1+theta^2*(nu^2-1))
190       
191        retval = 9*((sin(arg)-arg*cos(arg))/arg^3)^2
192       
193    return retval
194   
195End     //Function eor()
196
197// this is all there is to the smeared calculation!
198Function SmearedEllipsoidForm(s) :FitFunc
199        Struct ResSmearAAOStruct &s
200       
201        //the name of your unsmeared model is the first argument
202        Smear_Model_20(EllipsoidForm,s.coefW,s.xW,s.yW,s.resW)
203
204        return(0)
205End
206
207//wrapper to calculate the smeared model as an AAO-Struct
208// fills the struct and calls the ususal function with the STRUCT parameter
209//
210// used only for the dependency, not for fitting
211//
212Function fSmearedEllipsoidForm(coefW,yW,xW)
213        Wave coefW,yW,xW
214       
215        String str = getWavesDataFolder(yW,0)
216        String DF="root:"+str+":"
217       
218        WAVE resW = $(DF+str+"_res")
219       
220        STRUCT ResSmearAAOStruct fs
221        WAVE fs.coefW = coefW   
222        WAVE fs.yW = yW
223        WAVE fs.xW = xW
224        WAVE fs.resW = resW
225       
226        Variable err
227        err = SmearedEllipsoidForm(fs)
228       
229        return (0)
230End
Note: See TracBrowser for help on using the repository browser.