source: sans/Dev/trunk/NCNR_User_Procedures/SANS_Analysis/Models/UniformEllipsoid_v40.ipf @ 325

Last change on this file since 325 was 325, checked in by ajj, 15 years ago

Adding SANS Analysis Models to new Dev tree

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