source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/GaussSpheres.ipf @ 151

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

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

File size: 7.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "sphere"
5// plots the form factor of  spheres with a Gaussian radius distribution
6//
7// also can plot the distribution itself, based on the current model parameters
8//
9// integration is currently done using 20-pt quadrature, but may benefit from
10//switching to an adaptive integration.
11//
12
13Proc PlotGaussPolySphere(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_pgs,ywave_pgs
20        xwave_pgs = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
21        Make/O/D coef_pgs = {0.01,60,0.2,1e-6,3e-6,0.001}
22        make/O/T parameters_pgs = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
23        Edit parameters_pgs,coef_pgs
24       
25        Variable/G root:g_pgs
26        g_pgs := GaussPolySphere(coef_pgs,ywave_pgs,xwave_pgs)
27        Display ywave_pgs vs xwave_pgs
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// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
35Proc PlotSmearedGaussPolySpheres(str)                                                           
36        String str
37        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
38       
39        // if any of the resolution waves are missing => abort
40        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
41                Abort
42        endif
43       
44        SetDataFolder $("root:"+str)
45       
46        // Setup parameter table for model function
47        Make/O/D smear_coef_pgs = {0.01,60,0.2,1e-6,3e-6,0.001}                                 
48        make/o/t smear_parameters_pgs = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}     
49        Edit smear_parameters_pgs,smear_coef_pgs                                       
50       
51        // output smeared intensity wave, dimensions are identical to experimental QSIG values
52        // make extra copy of experimental q-values for easy plotting
53        Duplicate/O $(str+"_q") smeared_pgs,smeared_qvals                               
54        SetScale d,0,0,"1/cm",smeared_pgs                                                       
55                                       
56        Variable/G gs_pgs=0
57        gs_pgs := fSmearedGaussPolySphere(smear_coef_pgs,smeared_pgs,smeared_qvals)     //this wrapper fills the STRUCT
58       
59        Display smeared_pgs vs smeared_qvals                                                                   
60        ModifyGraph log=1,marker=29,msize=2,mode=4
61        Label bottom "q (\\S-1\\M)"
62        Label left "Intensity (cm\\S-1\\M)"
63        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
64       
65        SetDataFolder root:
66End
67       
68
69
70
71//AAO version, uses XOP if available
72// simply calls the original single point calculation with
73// a wave assignment (this will behave nicely if given point ranges)
74Function GaussPolySphere(cw,yw,xw) : FitFunc
75        Wave cw,yw,xw
76       
77#if exists("GaussPolySphereX")
78        yw = GaussPolySphereX(cw,xw)
79#else
80        yw = fGaussPolySphere(cw,xw)
81#endif
82        return(0)
83End
84
85Function fGaussPolySphere(w,xx) : FitFunc
86        wave w
87        variable xx
88       
89        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho
90       
91        //set up the coefficient values
92        scale=w[0]
93        rad=w[1]
94        pd=w[2]
95        sig=pd*rad
96        rho=w[3]
97        rhos=w[4]
98        delrho=rho-rhos
99        bkg=w[5]
100       
101//temp set scale=1 and bkg=0 for quadrature calc
102        Make/O/D/N=4 sphere_temp
103        sphere_temp[0] = 1
104        sphere_temp[1] = rad            //changed in loop
105        sphere_temp[2] = delrho
106        sphere_temp[3] = 0
107       
108        //could use 5 pt quadrature to integrate over the size distribution, since it's a gaussian
109        //currently using 20 pts...
110        Variable va,vb,ii,zi,nord,yy,summ,inten
111        String weightStr,zStr
112       
113        //select number of gauss points by setting nord=20 or76 points
114        nord = 20
115//      nord = 76
116       
117        weightStr = "gauss"+num2str(nord)+"wt"
118        zStr = "gauss"+num2str(nord)+"z"
119       
120        if (WaveExists($weightStr) == 0) // wave reference is not valid,
121                Make/D/N=(nord) $weightStr,$zStr
122                Wave gauWt = $weightStr
123                Wave gauZ = $zStr               // wave references to pass
124                if(nord==20)
125                        Make20GaussPoints(gauWt,gauZ)
126                else
127                        Make76GaussPoints(gauWt,gauZ)
128                endif   
129        else
130                if(exists(weightStr) > 1)
131                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
132                endif
133                Wave gauWt = $weightStr
134                Wave gauZ = $zStr               // create the wave references
135        endif
136       
137        // end points of integration
138        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
139        // +/- 3 sigq catches 99.73% of distrubution
140        // change limits (and spacing of zi) at each evaluation based on R()
141        //integration from va to vb
142        va = -4*sig + rad
143        if (va<0)
144                va=0            //to avoid numerical error when  va<0 (-ve q-value)
145        endif
146        vb = 4*sig +rad
147       
148        summ = 0.0              // initialize integral
149        Make/O/N=1 tmp_yw,tmp_xw
150        tmp_xw[0] = xx
151        for(ii=0;ii<nord;ii+=1)
152                // calculate Gauss points on integration interval (r-value for evaluation)
153                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
154                sphere_temp[1] = zi
155                // calculate sphere scattering
156                SphereForm(sphere_temp,tmp_yw,tmp_xw)           // AAO calculation, one point
157                yy = gauWt[ii] *  Gauss_distr(sig,rad,zi) * tmp_yw[0]
158                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
159               
160                summ += yy              //add to the running total of the quadrature
161        endfor
162// calculate value of integral to return
163        inten = (vb-va)/2.0*summ
164       
165        //re-normalize by polydisperse sphere volume
166        inten /= (4*pi/3*rad^3)*(1+3*pd^2)
167       
168        inten *= scale
169        inten+=bkg
170       
171        Return(inten)
172End
173
174//wrapper to calculate the smeared model as an AAO-Struct
175// fills the struct and calls the ususal function with the STRUCT parameter
176//
177// used only for the dependency, not for fitting
178//
179Function fSmearedGaussPolySphere(coefW,yW,xW)
180        Wave coefW,yW,xW
181       
182        String str = getWavesDataFolder(yW,0)
183        String DF="root:"+str+":"
184       
185        WAVE resW = $(DF+str+"_res")
186       
187        STRUCT ResSmearAAOStruct fs
188        WAVE fs.coefW = coefW   
189        WAVE fs.yW = yW
190        WAVE fs.xW = xW
191        WAVE fs.resW = resW
192       
193        Variable err
194        err = SmearedGaussPolySphere(fs)
195       
196        return (0)
197End
198
199// this is all there is to the smeared calculation!
200Function SmearedGaussPolySphere(s) :FitFunc
201        Struct ResSmearAAOStruct &s
202
203//      the name of your unsmeared model (AAO) is the first argument
204        Smear_Model_20(GaussPolySphere,s.coefW,s.xW,s.yW,s.resW)
205
206        return(0)
207End
208
209
210
211//calculates number density given the coefficients of the Gaussian distribution
212// the scale factor is the volume fraction
213// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
214Macro NumberDensity_Gauss()
215
216        Variable nden,vpoly,Ravg,p,Rg,I0,Sv,phi
217       
218        if(WaveExists(coef_pgs)==0)
219                abort "You need to plot the model first to create the coefficient table"
220        Endif
221        Ravg = coef_pgs[1]  // mean radius (A)
222        p = coef_pgs[2]  // polydispersity
223        phi = coef_pgs[0]  // volume fraction
224        Print "mean radius (A) = ",Ravg
225        Print "polydispersity (sig/avg) = ",p
226        Print "volume fraction = ",phi
227       
228        //re-normalize by polydisperse sphere volume
229        vpoly = (4*pi/3*Ravg^3)*(1+3*p^2)
230        nden = phi/vpoly                //nden in 1/A^3
231       
232        Print "number density (A^-3) = ",nden
233       
234        Rg = Ravg*(0.6)^0.5*(1+28*p^2+210*p^4+420*p^6+105*p^8) / (1+15*p^2+45*p^4+15*p^6)
235        Print "Guinier Radius (A) = ",Rg
236        I0 = 1.0e8*(4./3.)*pi*phi*(coef_pgs[3]- coef_pgs[4])^2*Ravg^3*(1+15*p^2+45*p^4+15*p^6)/(1+3*p^2)
237        Print "Forward scattering cross-section (cm-1 sr-1) I(0)= ",I0
238        Sv=1.0e8*(3*phi/Ravg)*(1+p^2)/(1+3*p^2)
239        Print "Interfacial surface area per unit sample volume (cm-1) Sv= ",Sv 
240        End
241
242// plots the Gauss distribution based on the coefficient values
243// a static calculation, so re-run each time
244//
245Macro PlotGaussDistribution()
246
247        variable pd,avg,zz,maxr
248       
249        if(WaveExists(coef_pgs)==0)
250                abort "You need to plot the model first to create the coefficient table"
251        Endif
252        pd=coef_pgs[2]
253        avg = coef_pgs[1]
254       
255        Make/O/D/N=1000 Gauss_distribution
256        maxr =  avg*(1+10*pd)
257
258        SetScale/I x, 0, maxr, Gauss_distribution
259        Gauss_distribution = Gauss_distr(pd*avg,avg,x)
260        Display Gauss_distribution
261        legend
262End
263
264
265Function Gauss_distr(sig,avg,pt)
266        Variable sig,avg,pt
267       
268        Variable retval
269       
270        retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2)
271        return(retval)
272End
Note: See TracBrowser for help on using the repository browser.