source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/SchulzSpheres_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: 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       
31        AddModelToStrings("SchulzSpheres","coef_sch","parameters_sch","sch")
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       
66        AddModelToStrings("SmearedSchulzSpheres","smear_coef_sch","smear_parameters_sch","sch")
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
182// add in the background
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//
297//      Variable scale,radius,pd,delrho,bkg,zz,rho,rhos
298//      scale = w[0]
299//      radius = w[1]
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
310//      Variable answer,zp1,zp2,zp3,vpoly
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
341//      a = radius*(1-range*pd)
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
370//      answer = (vb-va)/2.0*summ
371//     
372//      //re-normalize by the average volume
373//      zp1 = zz + 1.
374//              zp2 = zz + 2.
375//              zp3 = zz + 3.
376//              vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(radius)^3
377//      answer /= vpoly
378////scale
379//      answer *= scale
380//// add in the background
381//      answer += bkg
382//
383//      Return (answer)
384//End
385//
Note: See TracBrowser for help on using the repository browser.