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

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

Modified all procedure files to add to the keyword=value strings to identify the function, coefficients, and suffix once the model is plotted (and the objects will exist)

a one-liner in the Plot and PlotSmeared? macros.

  • necessary for smoother functioning of the wrapper panel.
File size: 7.9 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)
32       
33        AddModelToStrings("GaussSpheres","coef_pgs","pgs")
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 (\\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:
68        AddModelToStrings("SmearedGaussSpheres","smear_coef_pgs","pgs")
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       
92        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho
93       
94        //set up the coefficient values
95        scale=w[0]
96        rad=w[1]
97        pd=w[2]
98        sig=pd*rad
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=4 sphere_temp
106        sphere_temp[0] = 1
107        sphere_temp[1] = rad            //changed in loop
108        sphere_temp[2] = delrho
109        sphere_temp[3] = 0
110       
111        //could use 5 pt quadrature to integrate over the size distribution, since it's a gaussian
112        //currently using 20 pts...
113        Variable va,vb,ii,zi,nord,yy,summ,inten
114        String weightStr,zStr
115       
116        //select number of gauss points by setting nord=20 or76 points
117        nord = 20
118//      nord = 76
119       
120        weightStr = "gauss"+num2str(nord)+"wt"
121        zStr = "gauss"+num2str(nord)+"z"
122       
123        if (WaveExists($weightStr) == 0) // wave reference is not valid,
124                Make/D/N=(nord) $weightStr,$zStr
125                Wave gauWt = $weightStr
126                Wave gauZ = $zStr               // wave references to pass
127                if(nord==20)
128                        Make20GaussPoints(gauWt,gauZ)
129                else
130                        Make76GaussPoints(gauWt,gauZ)
131                endif   
132        else
133                if(exists(weightStr) > 1)
134                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
135                endif
136                Wave gauWt = $weightStr
137                Wave gauZ = $zStr               // create the wave references
138        endif
139       
140        // end points of integration
141        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
142        // +/- 3 sigq catches 99.73% of distrubution
143        // change limits (and spacing of zi) at each evaluation based on R()
144        //integration from va to vb
145        va = -4*sig + rad
146        if (va<0)
147                va=0            //to avoid numerical error when  va<0 (-ve q-value)
148        endif
149        vb = 4*sig +rad
150       
151        summ = 0.0              // initialize integral
152        Make/O/N=1 tmp_yw,tmp_xw
153        tmp_xw[0] = xx
154        for(ii=0;ii<nord;ii+=1)
155                // calculate Gauss points on integration interval (r-value for evaluation)
156                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
157                sphere_temp[1] = zi
158                // calculate sphere scattering
159                SphereForm(sphere_temp,tmp_yw,tmp_xw)           // AAO calculation, one point
160                yy = gauWt[ii] *  Gauss_distr(sig,rad,zi) * tmp_yw[0]
161                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
162               
163                summ += yy              //add to the running total of the quadrature
164        endfor
165// calculate value of integral to return
166        inten = (vb-va)/2.0*summ
167       
168        //re-normalize by polydisperse sphere volume
169        inten /= (4*pi/3*rad^3)*(1+3*pd^2)
170       
171        inten *= scale
172        inten+=bkg
173       
174        Return(inten)
175End
176
177//wrapper to calculate the smeared model as an AAO-Struct
178// fills the struct and calls the ususal function with the STRUCT parameter
179//
180// used only for the dependency, not for fitting
181//
182Function fSmearedGaussSpheres(coefW,yW,xW)
183        Wave coefW,yW,xW
184       
185        String str = getWavesDataFolder(yW,0)
186        String DF="root:"+str+":"
187       
188        WAVE resW = $(DF+str+"_res")
189       
190        STRUCT ResSmearAAOStruct fs
191        WAVE fs.coefW = coefW   
192        WAVE fs.yW = yW
193        WAVE fs.xW = xW
194        WAVE fs.resW = resW
195       
196        Variable err
197        err = SmearedGaussSpheres(fs)
198       
199        return (0)
200End
201
202// this is all there is to the smeared calculation!
203Function SmearedGaussSpheres(s) :FitFunc
204        Struct ResSmearAAOStruct &s
205
206//      the name of your unsmeared model (AAO) is the first argument
207        Smear_Model_20(GaussSpheres,s.coefW,s.xW,s.yW,s.resW)
208
209        return(0)
210End
211
212
213
214//calculates number density given the coefficients of the Gaussian distribution
215// the scale factor is the volume fraction
216// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
217Macro NumberDensity_Gauss()
218
219        Variable nden,vpoly,Ravg,p,Rg,I0,Sv,phi
220       
221        if(WaveExists(coef_pgs)==0)
222                abort "You need to plot the model first to create the coefficient table"
223        Endif
224        Ravg = coef_pgs[1]  // mean radius (A)
225        p = coef_pgs[2]  // polydispersity
226        phi = coef_pgs[0]  // volume fraction
227        Print "mean radius (A) = ",Ravg
228        Print "polydispersity (sig/avg) = ",p
229        Print "volume fraction = ",phi
230       
231        //re-normalize by polydisperse sphere volume
232        vpoly = (4*pi/3*Ravg^3)*(1+3*p^2)
233        nden = phi/vpoly                //nden in 1/A^3
234       
235        Print "number density (A^-3) = ",nden
236       
237        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)
238        Print "Guinier Radius (A) = ",Rg
239        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)
240        Print "Forward scattering cross-section (cm-1 sr-1) I(0)= ",I0
241        Sv=1.0e8*(3*phi/Ravg)*(1+p^2)/(1+3*p^2)
242        Print "Interfacial surface area per unit sample volume (cm-1) Sv= ",Sv 
243        End
244
245// plots the Gauss distribution based on the coefficient values
246// a static calculation, so re-run each time
247//
248Macro PlotGaussDistribution()
249
250        variable pd,avg,zz,maxr
251       
252        if(WaveExists(coef_pgs)==0)
253                abort "You need to plot the model first to create the coefficient table"
254        Endif
255        pd=coef_pgs[2]
256        avg = coef_pgs[1]
257       
258        Make/O/D/N=1000 Gauss_distribution
259        maxr =  avg*(1+10*pd)
260
261        SetScale/I x, 0, maxr, Gauss_distribution
262        Gauss_distribution = Gauss_distr(pd*avg,avg,x)
263        Display Gauss_distribution
264        legend
265End
266
267
268Function Gauss_distr(sig,avg,pt)
269        Variable sig,avg,pt
270       
271        Variable retval
272       
273        retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2)
274        return(retval)
275End
Note: See TracBrowser for help on using the repository browser.