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

Last change on this file since 999 was 999, checked in by srkline, 7 years ago

changes to a few analysis models to make these Igor 7-ready

adding mask editing utilities

many changes to event mode for easier processing of split lists

updated event mode help file

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