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

Last change on this file since 145 was 145, checked in by ajj, 15 years ago

Changes:

  • Made USANS weighting matrix double precision
  • Changed Smear_Model_(5,10,20,76) to be AAO and added test for presence of USANS matrix. If the number of columns in the resolution wave is > 4 then do the matrix multiplication, otherwise do gaussian quadrature via Smear_Model_N.
  • Changed Sphere.ipf and SchulzSpheres?.ipf to use AAO form of Smear_Model_20
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 PlotSchulzPolySpheres(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 PlotSmearedSchulzPolySpheres(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.