source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Fractal_v40.ipf @ 633

Last change on this file since 633 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.
File size: 5.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4// plots scattering from a mass fractal object
5// uses the model of Teixeria
6//
7// REFERENCE: J. Appl. Cryst. vol 21, p781-785
8//  Uses eq.1, 4, and 16
9//
10// Macro for fractal parameters added JGB 2004
11
12Proc PlotFractal(num,qmin,qmax)                                         
13        Variable num=128,qmin=0.001,qmax=0.5
14        Prompt num "Enter number of data points for model: "
15        Prompt qmin "Enter minimum q-value (A^-1) for model: "
16        Prompt qmax "Enter maximum q-value (A^-1) for model: "
17       
18        Make/O/D/n=(num) xwave_fra,ywave_fra                                   
19        xwave_fra = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))                                     
20        Make/O/D coef_fra = {0.05,5,2,100,2e-6,6.35e-6,0}                                               
21        make/o/t parameters_fra = {"Volume Fraction (scale)","Block Radius (A)","fractal dimension","correlation length (A)","SLD block (A-2)","SLD solvent (A-2)","bkgd (cm^-1 sr^-1)"}               
22        Edit parameters_fra,coef_fra           
23       
24        Variable/G root:g_fra
25        g_fra := Fractal(coef_fra,ywave_fra,xwave_fra)                 
26        Display ywave_fra vs xwave_fra                                                 
27        ModifyGraph log=1,marker=29,msize=2,mode=4                     
28        Label bottom "q (A\\S-1\\M)"
29        Label left "Intensity (cm\\S-1\\M)"                                     
30        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
31       
32        AddModelToStrings("Fractal","coef_fra","parameters_fra","fra")
33End
34
35// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
36Proc PlotSmearedFractal(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_fra =  {0.05,5,2,100,2e-6,6.35e-6,0}
49        make/o/t smear_parameters_fra = {"Volume Fraction (scale)","Block Radius (A)","fractal dimension","correlation length (A)","SLD block (A-2)","SLD solvent (A-2)","bkgd (cm^-1 sr^-1)"}
50        Edit smear_parameters_fra,smear_coef_fra                                        //display parameters in a table
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_fra,smeared_qvals                               //
55        SetScale d,0,0,"1/cm",smeared_fra                                               //
56                                       
57        Variable/G gs_fra=0
58        gs_fra := fSmearedFractal(smear_coef_fra,smeared_fra,smeared_qvals)     //this wrapper fills the STRUCT
59       
60        Display smeared_fra vs smeared_qvals           
61        ModifyGraph log=1,marker=29,msize=2,mode=4
62        Label bottom "q (A\\S-1\\M)"
63        Label left "I(q) (cm\\S-1\\M)"
64        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
65       
66        SetDataFolder root:
67        AddModelToStrings("SmearedFractal","smear_coef_fra","smear_parameters_fra","fra")
68End
69
70
71
72//calculates the physical parameters related to the
73//model parameters. See the reference at the top of the
74//file for details
75Macro NumberDensity_Fractal()
76       
77        Variable nden,phi,r0,Df,corr,s0,vpoly,i0,rg
78       
79        if(Exists("coef_fra")!=1)
80                abort "You need to plot the unsmeared model first to create the coefficient table"
81        Endif
82       
83        phi = coef_fra[0]   // volume fraction of building blocks
84        r0 = coef_fra[1]    // building block radius
85        Df = coef_fra[2]    // fractal dimension
86        corr = coef_fra[3]  // fractal correlation length (of cluster)
87       
88        Print "mean building block radius (A) = ",r0
89        Print "volume fraction = ",phi
90       
91 //  average particle volume   
92        vpoly = 4*Pi/3*r0^3
93        nden = phi/vpoly                //nden in 1/A^3
94        i0 = 1.0e8*phi*vpoly*(coef_fra[4]-coef_fra[5])^2  // 1/cm/sr
95        rg = corr*( Df*(Df+1)/2 )^0.5
96        s0 = exp(gammln(Df+1))*(corr/r0)^Df
97        Print "number density (A^-3) = ",nden
98        Print "Guinier radius (A) = ",rg
99        Print "Aggregation number G = ",s0
100        Print "Forward cross section of building blocks (cm-1 sr-1) I(0) = ",i0
101        Print "Forward cross section of clusters (cm-1 sr-1) I(0) = ",i0*s0
102End
103
104//AAO version, uses XOP if available
105// simply calls the original single point calculation with
106// a wave assignment (this will behave nicely if given point ranges)
107Function Fractal(cw,yw,xw) : FitFunc
108        Wave cw,yw,xw
109       
110#if exists("FractalX")
111        yw = FractalX(cw,xw)
112#else
113        yw = fFractal(cw,xw)
114#endif
115        return(0)
116End
117
118//fractal scattering function
119Function fFractal(w,x) :FitFunc
120        wave w
121        variable x
122       
123        variable r0,Df,corr,phi,sldp,sldm,bkg,bes
124        variable pq,sq,ans
125        phi=w[0]   // volume fraction of building block spheres...
126        r0=w[1]    //  radius of building block
127        Df=w[2]     //  fractal dimension
128        corr=w[3] //  correlation length of fractal-like aggregates
129        sldp = w[4] // SLD of building block
130        sldm = w[5] // SLD of matrix or solution
131        bkg=w[6]  //  flat background
132       
133        //calculate P(q) for the spherical subunits, units cm-1 sr-1
134        if(x*r0 == 0)
135                bes = 1
136        else
137                bes = (3*(sin(x*r0) - x*r0*cos(x*r0))/(x*r0)^3)^2
138        endif
139       
140        pq = 1.0e8*phi*4/3*pi*r0^3*(sldp-sldm)^2*bes
141
142        //calculate S(q)
143        sq = Df*exp(gammln(Df-1))*sin((Df-1)*atan(x*corr))
144        sq /= (x*r0)^Df * (1 + 1/(x*corr)^2)^((Df-1)/2)
145        sq += 1
146
147        //combine and return
148        ans = pq*sq + bkg
149
150        return (ans)
151End
152
153//wrapper to calculate the smeared model as an AAO-Struct
154// fills the struct and calls the ususal function with the STRUCT parameter
155//
156// used only for the dependency, not for fitting
157//
158Function fSmearedFractal(coefW,yW,xW)
159        Wave coefW,yW,xW
160       
161        String str = getWavesDataFolder(yW,0)
162        String DF="root:"+str+":"
163       
164        WAVE resW = $(DF+str+"_res")
165       
166        STRUCT ResSmearAAOStruct fs
167        WAVE fs.coefW = coefW   
168        WAVE fs.yW = yW
169        WAVE fs.xW = xW
170        WAVE fs.resW = resW
171       
172        Variable err
173        err = SmearedFractal(fs)
174       
175        return (0)
176End
177
178//the smeared model calculation
179Function SmearedFractal(s) :FitFunc
180        Struct ResSmearAAOStruct &s
181
182//      the name of your unsmeared model (AAO) is the first argument
183        Smear_Model_20(Fractal,s.coefW,s.xW,s.yW,s.resW)
184
185        return(0)
186End
Note: See TracBrowser for help on using the repository browser.