source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/GaussSpheres.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: 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 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 (^-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 := 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 (\\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 PlotSmearedGaussSpheres(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 := fSmearedGaussSpheres(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 GaussSpheres(cw,yw,xw) : FitFunc
75        Wave cw,yw,xw
76       
77#if exists("GaussSpheresX")
78        yw = GaussSpheresX(cw,xw)
79#else
80        yw = fGaussSpheres(cw,xw)
81#endif
82        return(0)
83End
84
85Function fGaussSpheres(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 fSmearedGaussSpheres(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 = SmearedGaussSpheres(fs)
195       
196        return (0)
197End
198
199// this is all there is to the smeared calculation!
200Function SmearedGaussSpheres(s) :FitFunc
201        Struct ResSmearAAOStruct &s
202
203//      the name of your unsmeared model (AAO) is the first argument
204        Smear_Model_20(GaussSpheres,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.