source:sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/SchulzSpheres_v40.ipf@379

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

File size: 10.2 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 for spheres with a Sculz distribution of radius
6// at low polydispersity (< 0.2), it is very similar to the Gaussian distribution
7// at larger polydispersities, it is more skewed and similar to log-normal
8//
9
10//
11Proc PlotSchulzSpheres(num,qmin,qmax)
12        Variable num=128,qmin=0.001,qmax=0.7
13        Prompt num "Enter number of data points for model: "
14        Prompt qmin "Enter minimum q-value (A^-1) for model: "
15        Prompt qmax "Enter maximum q-value (A^-1) for model: "
16
17        Make/O/D/N=(num) xwave_sch,ywave_sch
18        xwave_sch = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
19        Make/O/D coef_sch = {0.01,60,0.2,1e-6,3e-6,0.001}
20        make/O/T parameters_sch = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
21        Edit parameters_sch,coef_sch
22
23        Variable/G root:g_sch
24        g_sch := SchulzSpheres(coef_sch,ywave_sch,xwave_sch)
25        Display ywave_sch vs xwave_sch
26        ModifyGraph log=1,marker=29,msize=2,mode=4
27        Label bottom "q (A\\S-1\\M)"
28        Label left "Intensity (cm\\S-1\\M)"
29        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
30
32End
33
34
35// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
36Proc PlotSmearedSchulzSpheres(str)
37        String str
38        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
39
40        // if any of the resolution waves are missing => abort
41        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
42                Abort
43        endif
44
45        SetDataFolder \$("root:"+str)
46
47        // Setup parameter table for model function
48        Make/O/D smear_coef_sch = {0.01,60,0.2,1e-6,3e-6,0.001}
49        make/o/t smear_parameters_sch = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
50        Edit smear_parameters_sch,smear_coef_sch
51
52        // output smeared intensity wave, dimensions are identical to experimental QSIG values
53        // make extra copy of experimental q-values for easy plotting
54        Duplicate/O \$(str+"_q") smeared_sch,smeared_qvals
55        SetScale d,0,0,"1/cm",smeared_sch
56
57        Variable/G gs_sch=0
58        gs_sch := fSmearedSchulzSpheres(smear_coef_sch,smeared_sch,smeared_qvals)       //this wrapper fills the STRUCT
59
60        Display smeared_sch vs smeared_qvals
61        ModifyGraph log=1,marker=29,msize=2,mode=4
62        Label bottom "q (A\\S-1\\M)"
63        Label left "Intensity (cm\\S-1\\M)"
64        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
65
67End
68
69
70
71Function Schulz_Point(x,avg,zz)
72        Variable x,avg,zz
73
74        Variable dr
75
76        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
77        return (exp(dr))
78End
79
80
81//AAO version, uses XOP if available
82// simply calls the original single point calculation with
83// a wave assignment (this will behave nicely if given point ranges)
84Function SchulzSpheres(cw,yw,xw) : FitFunc
85        Wave cw,yw,xw
86
87#if exists("SchulzSpheresX")
88        yw = SchulzSpheresX(cw,xw)
89#else
90        yw = fSchulzSpheres(cw,xw)
91#endif
92        return(0)
93End
94
95//use the analytic formula from Kotlarchyk & Chen, JCP 79 (1983) 2461
96//equations 23-30
97//
98// need to calculate in terms of logarithms to avoid numerical errors
99//
100Function fSchulzSpheres(w,x) : FitFunc
101        Wave w
102        Variable x
103
104        Variable scale,ravg,pd,delrho,bkg,zz,rho,rhos,vpoly
105        scale = w[0]
106        ravg = w[1]
107        pd = w[2]
108        rho = w[3]
109        rhos = w[4]
110        bkg = w[5]
111
112        delrho=rho-rhos
113        zz = (1/pd)^2-1
114
115        Variable zp1,zp2,zp3,zp4,zp5,zp6,zp7
116        Variable aa,b1,b2,b3,at1,at2,rt1,rt2,rt3,t1,t2,t3
117        Variable v1,v2,v3,g1,g11,gd,pq,g2,g22
118
119        ZP1 = zz + 1
120        ZP2 = zz + 2
121        ZP3 = zz + 3
122        ZP4 = zz + 4
123        ZP5 = zz + 5
124        ZP6 = zz + 6
125        ZP7 = zz + 7
126
127        //small QR limit - use Guinier approx
128        Variable i_zero,Rg2,zp8
129        zp8 = zz+8
130        if(x*ravg < 0.1)
131                i_zero = scale*delrho*delrho*1e8*4*Pi/3*ravg^3
132                i_zero *= zp6*zp5*zp4/zp1/zp1/zp1               //6th moment / 3rd moment
133                Rg2 = 3*zp8*zp7/5/(zp1^2)*ravg*ravg
134                pq = i_zero*exp(-x*x*Rg2/3)
135                pq += bkg
136                return(pq)
137        endif
138//
139        aa = (zz+1)/x/Ravg
140
141        AT1 = atan(1/aa)
142        AT2 = atan(2/aa)
143//
144//  CALCULATIONS ARE PERFORMED TO AVOID  LARGE # ERRORS
145// - trick is to propogate the a^(z+7) term through the G1
146//
147        T1 = ZP7*log(aa) - zp1/2*log(aa*aa+4)
148        T2 = ZP7*log(aa) - zp3/2*log(aa*aa+4)
149        T3 = ZP7*log(aa) - zp2/2*log(aa*aa+4)
150//      Print T1,T2,T3
151        RT1 = alog(T1)
152        RT2 = alog(T2)
153        RT3 = alog(T3)
154        V1 = aa^6 - RT1*cos(zp1*at2)
155        V2 = ZP1*ZP2*( aa^4 + RT2*cos(zp3*at2) )
156        V3 = -2*ZP1*RT3*SIN(zp2*at2)
157        G1 = (V1+V2+V3)
158
159        Pq = log(G1) - 6*log(ZP1) + 6*log(Ravg)
160        Pq = alog(Pq)*8*PI*PI*delrho*delrho
161
162//
163// beta factor is not used here, but could be for the
164// decoupling approximation
165//
166//      G11 = G1
167//      GD = -ZP7*log(aa)
168//      G1 = log(G11) + GD
169//
170//      T1 = ZP1*at1
171//      T2 = ZP2*at1
172//      G2 = SIN( T1 ) - ZP1/SQRT(aa*aa+1)*COS( T2 )
173//      G22 = G2*G2
174//      BETA = ZP1*log(aa) - ZP1*log(aa*aa+1) - G1 + log(G22)
175//      BETA = 2*alog(BETA)
176
177//re-normalize by the average volume
178        vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(ravg)^3
179        Pq /= vpoly
180//scale, convert to cm^-1
181        Pq *= scale * 1e8
183        Pq += bkg
184
185        //return (g1)
186        Return (Pq)
187End
188
189//wrapper to calculate the smeared model as an AAO-Struct
190// fills the struct and calls the ususal function with the STRUCT parameter
191//
192// used only for the dependency, not for fitting
193//
194Function fSmearedSchulzSpheres(coefW,yW,xW)
195        Wave coefW,yW,xW
196
197        String str = getWavesDataFolder(yW,0)
198        String DF="root:"+str+":"
199
200        WAVE resW = \$(DF+str+"_res")
201
202        STRUCT ResSmearAAOStruct fs
203        WAVE fs.coefW = coefW
204        WAVE fs.yW = yW
205        WAVE fs.xW = xW
206        WAVE fs.resW = resW
207
208        Variable err
209        err = SmearedSchulzSpheres(fs)
210
211        return (0)
212End
213
214// this is all there is to the smeared calculation!
215Function SmearedSchulzSpheres(s) :FitFunc
216        Struct ResSmearAAOStruct &s
217
218//      the name of your unsmeared model (AAO) is the first argument
219        Smear_Model_20(SchulzSpheres,s.coefW,s.xW,s.yW,s.resW)
220
221        return(0)
222End
223
224
225//calculates number density given the coefficients of the Schulz distribution
226// the scale factor is the volume fraction
227// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
228Macro NumberDensity_Schulz()
229
230        Variable nden,zz,zp1,zp2,zp3,zp4,zp5,zp6,zp7,zp8,vpoly,v2poly,rg,sv,i0
231
232        if(Exists("coef_sch")!=1)
233                abort "You need to plot the unsmeared model first to create the coefficient table"
234        Endif
235
236        Print "mean radius (A) = ",coef_sch[1]
237        Print "polydispersity (sig/avg) = ",coef_sch[2]
238        Print "volume fraction = ",coef_sch[0]
239
240        zz = (1/coef_sch[2])^2-1
241        zp1 = zz + 1.
242        zp2 = zz + 2.
243        zp3 = zz + 3.
244        zp4 = zz + 4.
245        zp5 = zz + 5.
246        zp6 = zz + 6.
247        zp7 = zz + 7.
248        zp8 = zz + 8.
249 //  average particle volume
250        vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(coef_sch[1])^3
251 //  average particle volume
252        v2poly = (4*Pi/3)^2*zp6*zp5*zp4*zp3*zp2/(zp1^5)*(coef_sch[1])^6
253        nden = coef_sch[0]/vpoly                //nden in 1/A^3
254        rg = coef_sch[1]*((3*zp8*zp7)/5/zp1/zp1)^0.5   // in A
255        sv = 1.0e8*3*coef_sch[0]*zp1/coef_sch[1]/zp3  // in 1/cm
256        i0 = 1.0e8*nden*v2poly*(coef_sch[3]-coef_sch[4])^2  // 1/cm/sr
257
258        Print "number density (A^-3) = ",nden
259        Print "Guinier radius (A) = ",rg
260        Print "Interfacial surface area / volume (cm-1) Sv = ",sv
261        Print "Forward cross section (cm-1 sr-1) I(0) = ",i0
262End
263
264// plots the Schulz distribution based on the coefficient values
265// a static calculation, so re-run each time
266//
267Macro PlotSchulzDistribution()
268
269        variable pd,avg,zz,maxr
270
271        if(Exists("coef_sch")!=1)
272                abort "You need to plot the unsmeared model first to create the coefficient table"
273        Endif
274        pd=coef_sch[2]
275        avg = coef_sch[1]
276        zz = (1/pd)^2-1
277
278        Make/O/D/N=1000 Schulz_distribution
279        maxr =  avg*(1+10*pd)
280
281        SetScale/I x, 0, maxr, Schulz_distribution
282        Schulz_distribution = Schulz_Point(x,avg,zz)
283        Display Schulz_distribution
284        legend
285End
286
287// don't use the integral technique here since there's an analytic solution
288// available
289//////////////////////
290//requires that sphere.ipf is included
291//
292//
293//Function SchulzSpheres_Integrated(w,x)
294//      Wave w
295//      Variable x
296//
298//      scale = w[0]
300//      pd = w[2]
301//      rho = w[3]
302//      rhos = w[4]
303//      bkg = w[5]
304//
305//      delrho=rho-rhos
306//      zz = (1/pd)^2-1
307////
308//// local variables
309//      Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
311//      String weightStr,zStr
312//
313//      //select number of gauss points by setting nord=20 or76 points
314////    nord = 20
315//      nord = 76
316//
317//      weightStr = "gauss"+num2str(nord)+"wt"
318//      zStr = "gauss"+num2str(nord)+"z"
319//
320//      if (WaveExists(\$weightStr) == 0) // wave reference is not valid,
321//              Make/D/N=(nord) \$weightStr,\$zStr
322//              Wave gauWt = \$weightStr
323//              Wave gauZ = \$zStr               // wave references to pass
324//              if(nord==20)
325//                      Make20GaussPoints(gauWt,gauZ)
326//              else
327//                      Make76GaussPoints(gauWt,gauZ)
328//              endif
329//      else
330//              if(exists(weightStr) > 1)
331//                       Abort "wave name is already in use"            //executed only if name is in use elsewhere
332//              endif
333//              Wave gauWt = \$weightStr
334//              Wave gauZ = \$zStr               // create the wave references
335//      endif
336//
337//// set up the integration
338//// end points and weights
339//// limits are technically 0-inf, but wisely choose non-zero region of distribution
340//      Variable range=8                //multiples of the std. dev. from the mean
342//      if (a<0)
343//              a=0             //otherwise numerical error when pd >= 0.3, making a<0
344//      endif
345//      If(pd>0.3)
346//              range = 3.4 + (pd-0.3)*18               //to account for skewed tail
347//      Endif
348//      b = radius*(1+range*pd) // is this far enough past avg radius?
349//      va =a
350//      vb =b
351//
352////temp set scale=1 and bkg=0 for quadrature calc
353//      Make/O/D/N=4 sphere_temp
354//      sphere_temp[0] = 1
355//      sphere_temp[1] = radius         //changed in loop
356//      sphere_temp[2] = delrho
357//      sphere_temp[3] = 0
358//
359//// evaluate at Gauss points
360//      summ = 0.0              // initialize integral
361//      for(ii=0;ii<nord;ii+=1)
362//              zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
363//              sphere_temp[1] = zi
364//              yyy = gauWt[ii] * Schulz_Point(zi,radius,zz) * SphereForm(sphere_temp,x)
365//              //un-normalize by volume
366//              yyy *= 4*pi/3*zi^3
367//              summ += yyy
368//      endfor
369//// calculate value of integral to return
371//
372//      //re-normalize by the average volume
373//      zp1 = zz + 1.
374//              zp2 = zz + 2.
375//              zp3 = zz + 3.