# source:sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/GaussSpheres_v40.ipf@381

Last change on this file since 381 was 381, checked in by srkline, 14 years ago

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

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