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

Last change on this file since 835 was 835, checked in by srkline, 11 years ago

moving the Polarization routines to a temporary "beta" menu location on the macros menu, numbering the panels.

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