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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

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