source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/GaussSpheres_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 8.0 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
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 (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       
33        AddModelToStrings("GaussSpheres","coef_pgs","parameters_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 (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:
68        AddModelToStrings("SmearedGaussSpheres","smear_coef_pgs","smear_parameters_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=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
150        vb = 4*sig +rad
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
170        inten /= (4*pi/3*rad^3)*(1+3*pd^2)
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.