source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/LogNormalSphere_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.5 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 spherical particles with a log-normal radius distribution
6//
7// for the integration it may be better to use adaptive routine
8
9//Proc to setup data and coefficients to plot the model
10Proc PlotLogNormalSphere(num,qmin,qmax)
11        Variable num=128,qmin=0.001,qmax=0.7
12        Prompt num "Enter number of data points for model: "
13        Prompt qmin "Enter minimum q-value (A^-1) for model: "
14        Prompt qmax "Enter maximum q-value (A^-1) for model: "
15       
16        Make/O/D/N=(num) xwave_lns,ywave_lns
17        xwave_lns = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
18        Make/O/D coef_lns = {0.01,60,0.2,1e-6,2e-6,0}
19        make/O/T parameters_lns = {"Volume Fraction (scale)","exp(mu)=median Radius (A)","sigma","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
20        Edit parameters_lns,coef_lns
21       
22        Variable/G root:g_lns
23        g_lns := LogNormalSphere(coef_lns,ywave_lns,xwave_lns)
24        Display ywave_lns vs xwave_lns
25        ModifyGraph log=1,marker=29,msize=2,mode=4
26        Label bottom "q (A\\S-1\\M)"
27        Label left "Intensity (cm\\S-1\\M)"
28        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
29       
30        AddModelToStrings("LogNormalSphere","coef_lns","parameters_lns","lns")
31End
32
33// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
34Proc PlotSmearedLogNormalSphere(str)                                                           
35        String str
36        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
37       
38        // if any of the resolution waves are missing => abort
39        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
40                Abort
41        endif
42       
43        SetDataFolder $("root:"+str)
44       
45        // Setup parameter table for model function
46        Make/O/D smear_coef_lns = {0.01,60,0.2,1e-6,2e-6,0}                                     
47        make/o/t smear_parameters_lns = {"Volume Fraction (scale)","exp(mu)=median Radius (A)","sigma","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}               
48        Edit smear_parameters_lns,smear_coef_lns                                       
49       
50        // output smeared intensity wave, dimensions are identical to experimental QSIG values
51        // make extra copy of experimental q-values for easy plotting
52        Duplicate/O $(str+"_q") smeared_lns,smeared_qvals                               
53        SetScale d,0,0,"1/cm",smeared_lns                                                       
54                                       
55        Variable/G gs_lns=0
56        gs_lns := fSmearedLogNormalSphere(smear_coef_lns,smeared_lns,smeared_qvals)     //this wrapper fills the STRUCT
57       
58        Display smeared_lns vs smeared_qvals                                                                   
59        ModifyGraph log=1,marker=29,msize=2,mode=4
60        Label bottom "q (A\\S-1\\M)"
61        Label left "Intensity (cm\\S-1\\M)"
62        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
63       
64        SetDataFolder root:
65        AddModelToStrings("SmearedLogNormalSphere","smear_coef_lns","smear_parameters_lns","lns")
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 LogNormalSphere(cw,yw,xw) : FitFunc
75        Wave cw,yw,xw
76       
77#if exists("LogNormalSphereX")
78        yw = LogNormalSphereX(cw,xw)
79#else
80        yw = fLogNormalSphere(cw,xw)
81#endif
82        return(0)
83End
84
85
86// calculates the model at each q-value by integrating over the normalized size distribution
87// integration is done by gauss quadrature of either 20 or 76 points (nord)
88// 76 points is slower, but reccommended to remove high-q oscillations
89//
90Function fLogNormalSphere(w,xx): FitFunc
91        wave w
92        variable xx
93       
94        Variable scale,rad,sig,rho,rhos,bkg,delrho,mu,r3
95       
96        //set up the coefficient values
97        scale=w[0]
98        rad=w[1]                //rad is the median radius
99        mu = ln(w[1])
100        sig=w[2]
101        rho=w[3]
102        rhos=w[4]
103        delrho=rho-rhos
104        bkg=w[5]
105       
106//temp set scale=1 and bkg=0 for quadrature calc
107        Make/O/D/N=4 sphere_temp
108        sphere_temp[0] = 1
109        sphere_temp[1] = rad            //changed in loop
110        sphere_temp[2] = rho
111        sphere_temp[3] = rhos
112        sphere_temp[4] = 0
113       
114        //currently using 20 pts...
115        Variable va,vb,ii,zi,nord,yy,summ,inten
116        String weightStr,zStr
117       
118        //select number of gauss points by setting nord=20 or76 points
119//      nord = 20
120        nord = 76
121       
122        weightStr = "gauss"+num2str(nord)+"wt"
123        zStr = "gauss"+num2str(nord)+"z"
124       
125        if (WaveExists($weightStr) == 0) // wave reference is not valid,
126                Make/D/N=(nord) $weightStr,$zStr
127                Wave gauWt = $weightStr
128                Wave gauZ = $zStr               // wave references to pass
129                if(nord==20)
130                        Make20GaussPoints(gauWt,gauZ)
131                else
132                        Make76GaussPoints(gauWt,gauZ)
133                endif   
134        else
135                if(exists(weightStr) > 1)
136                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
137                endif
138                Wave gauWt = $weightStr
139                Wave gauZ = $zStr               // create the wave references
140        endif
141       
142        // end points of integration
143        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
144        // +/- 3 sigq catches 99.73% of distrubution
145        // change limits (and spacing of zi) at each evaluation based on R()
146        //integration from va to vb
147       
148//      va = -3*sig + rad
149        va = -3.5*sig +mu               //in ln(R) space
150        va = exp(va)                    //in R space
151        if (va<0)
152                va=0            //to avoid numerical error when  va<0 (-ve q-value)
153        endif
154//      vb = 3*exp(sig) +rad
155        vb = 3.5*sig*(1+sig)+ mu
156        vb = exp(vb)
157       
158        summ = 0.0              // initialize integral
159        Make/O/N=1 tmp_yw,tmp_xw
160        tmp_xw[0] = xx
161        for(ii=0;ii<nord;ii+=1)
162                // calculate Gauss points on integration interval (r-value for evaluation)
163                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
164                sphere_temp[1] = zi
165                // calculate sphere scattering
166                SphereForm(sphere_temp,tmp_yw,tmp_xw)                   //AAO calculation, one point wave
167                yy = gauWt[ii] *  LogNormal_distr(sig,mu,zi) * tmp_yw[0]
168                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
169               
170                summ += yy              //add to the running total of the quadrature
171        endfor
172// calculate value of integral to return
173        inten = (vb-va)/2.0*summ
174       
175        //re-normalize by polydisperse sphere volume
176        //third moment
177        r3 = exp(3*mu + 9/2*sig^2)              // <R^3> directly
178        inten /= (4*pi/3*r3)            //polydisperse volume
179       
180        inten *= scale
181        inten+=bkg
182       
183        Return(inten)
184End
185
186//wrapper to calculate the smeared model as an AAO-Struct
187// fills the struct and calls the ususal function with the STRUCT parameter
188//
189// used only for the dependency, not for fitting
190//
191Function fSmearedLogNormalSphere(coefW,yW,xW)
192        Wave coefW,yW,xW
193       
194        String str = getWavesDataFolder(yW,0)
195        String DF="root:"+str+":"
196       
197        WAVE resW = $(DF+str+"_res")
198       
199        STRUCT ResSmearAAOStruct fs
200        WAVE fs.coefW = coefW   
201        WAVE fs.yW = yW
202        WAVE fs.xW = xW
203        WAVE fs.resW = resW
204       
205        Variable err
206        err = SmearedLogNormalSphere(fs)
207       
208        return (0)
209End
210
211// this is all there is to the smeared calculation!
212Function SmearedLogNormalSphere(s) :FitFunc
213        Struct ResSmearAAOStruct &s
214
215//      the name of your unsmeared model (AAO) is the first argument
216        Smear_Model_20(LogNormalSphere,s.coefW,s.xW,s.yW,s.resW)
217
218        return(0)
219End
220       
221
222
223// normalization is correct, using 3rd moment of lognormal distribution
224//
225Function LogNormal_distr(sig,mu,pt)
226        Variable sig,mu,pt
227       
228        Variable retval
229       
230        retval = (1/ ( sig*pt*sqrt(2*pi)) )*exp( -0.5*((ln(pt) - mu)^2)/sig^2 )
231        return(retval)
232End
233
234//calculates number density given the coefficients of the lognormal distribution
235// the scale factor is the volume fraction
236// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
237Macro NumberDensity_LogN()
238       
239        Variable nden,r3,rg,sv,i0,ravg,rpk
240       
241        if(Exists("coef_lns")!=1)
242                abort "You need to plot the unsmeared model first to create the coefficient table"
243        Endif
244       
245        Print "median radius (A) = ",coef_lns[1]
246        Print "sigma = ",coef_lns[2]
247        Print "volume fraction = ",coef_lns[0]
248       
249        r3 = exp(3*ln(coef_lns[1]) + 9/2*coef_lns[2]^2)         // <R^3> directly,[A^3]
250        nden = coef_lns[0]/(4*pi/3*r3)          //nden in 1/A^3
251        ravg = exp(ln(coef_lns[1]) + 0.5*coef_lns[2]^2)
252        rpk = exp(ln(coef_lns[1]) - coef_lns[2]^2)
253        rg = (3./5.)^0.5*exp(ln(coef_lns[1]) + 7.*coef_lns[2]^2)
254        sv = 1.0e8*3*coef_lns[0]*exp(-ln(coef_lns[1]) - 2.5*coef_lns[2]^2)
255        i0 = 1.0e8*(4*pi/3)*coef_lns[0]*(coef_lns[3]-coef_lns[4])^2*exp(3*ln(coef_lns[1]) + 13.5*coef_lns[2]^2)
256       
257        Print "number density (A^-3) = ",nden
258        Print "mean radius (A) = ",ravg
259        Print "peak dis. radius (A) = ",rpk
260        Print "Guinier radius (A) = ",rg
261        Print "Interfacial surface area / volume (cm-1) Sv = ",sv
262        Print "Forward cross section (cm-1 sr-1) I(0) = ",i0
263End
264
265// plots the lognormal distribution based on the coefficient values
266// a static calculation, so re-run each time
267//
268Macro PlotLogNormalDistribution()
269
270        variable sig,mu,maxr
271       
272        if(Exists("coef_lns")!=1)
273                abort "You need to plot the unsmeared model first to create the coefficient table"
274        Endif
275        sig=coef_lns[2]
276        mu = ln(coef_lns[1])
277       
278        Make/O/D/N=1000 lognormal_distribution
279        maxr = 5*sig*(1+sig)+ mu
280        maxr = exp(maxr)
281        SetScale/I x, 0, maxr, lognormal_distribution
282        lognormal_distribution = LogNormal_distr(sig,mu,x)
283        Display lognormal_distribution
284        modifygraph log(bottom)=1
285        legend
286End
Note: See TracBrowser for help on using the repository browser.