source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/Fractal_PolyCore_v40.ipf @ 570

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

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

File size: 6.8 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// - basic subunit is a polydisperse (Schulz) sphere w/ shell - SRK Jun 2009
11// - based on two existing XOPs, so no need to write a new one
12//
13// Macro for fractal parameters added JGB 2004
14
15Proc PlotFractalPolyCore(num,qmin,qmax)                                         
16        Variable num=128,qmin=0.001,qmax=0.5
17        Prompt num "Enter number of data points for model: "
18        Prompt qmin "Enter minimum q-value (A^-1) for model: "
19        Prompt qmax "Enter maximum q-value (A^-1) for model: "
20       
21        Make/O/D/n=(num) xwave_fraPolyCore,ywave_fraPolyCore                                   
22        xwave_fraPolyCore = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))                                     
23        Make/O/D coef_fraPolyCore = {0.05,20,0.1,5,2,100,3.5e-6,1e-6,6.35e-6,0}                                         
24        make/o/t parameters_fraPolyCore = {"Volume Fraction (scale)","Block Radius (A)","block polydispersity (0,1)","shell thickness (A)","fractal dimension","correlation length (A)","SLD block (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkgd (cm^-1 sr^-1)"}           
25        Edit parameters_fraPolyCore,coef_fraPolyCore           
26       
27        Variable/G root:g_fraPolyCore
28        g_fraPolyCore := FractalPolyCore(coef_fraPolyCore,ywave_fraPolyCore,xwave_fraPolyCore)                 
29        Display ywave_fraPolyCore vs xwave_fraPolyCore                                                 
30        ModifyGraph log=1,marker=29,msize=2,mode=4                     
31        Label bottom "q (A\\S-1\\M)"
32        Label left "Intensity (cm\\S-1\\M)"                                     
33        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
34       
35        AddModelToStrings("FractalPolyCore","coef_fraPolyCore","parameters_fraPolyCore","fraPolyCore")
36End
37
38// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
39Proc PlotSmearedFractalPolyCore(str)                                                           
40        String str
41        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
42       
43        // if any of the resolution waves are missing => abort
44        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
45                Abort
46        endif
47       
48        SetDataFolder $("root:"+str)
49       
50        // Setup parameter table for model function
51        Make/O/D smear_coef_fraPolyCore =  {0.05,20,0.1,5,2,100,3.5e-6,1e-6,6.35e-6,0}
52        make/o/t smear_parameters_fraPolyCore = {"Volume Fraction (scale)","Block Radius (A)","block polydispersity (0,1)","shell thickness (A)","fractal dimension","correlation length (A)","SLD block (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkgd (cm^-1 sr^-1)"}
53        Edit smear_parameters_fraPolyCore,smear_coef_fraPolyCore                                        //display parameters in a table
54       
55        // output smeared intensity wave, dimensions are identical to experimental QSIG values
56        // make extra copy of experimental q-values for easy plotting
57        Duplicate/O $(str+"_q") smeared_fraPolyCore,smeared_qvals                               //
58        SetScale d,0,0,"1/cm",smeared_fraPolyCore                                               //
59                                       
60        Variable/G gs_fraPolyCore=0
61        gs_fraPolyCore := fSmearedFractalPolyCore(smear_coef_fraPolyCore,smeared_fraPolyCore,smeared_qvals)     //this wrapper fills the STRUCT
62       
63        Display smeared_fraPolyCore vs smeared_qvals           
64        ModifyGraph log=1,marker=29,msize=2,mode=4
65        Label bottom "q (A\\S-1\\M)"
66        Label left "I(q) (cm\\S-1\\M)"
67        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
68       
69        SetDataFolder root:
70        AddModelToStrings("SmearedFractalPolyCore","smear_coef_fraPolyCore","smear_parameters_fraPolyCore","fraPolyCore")
71End
72
73
74
75//calculates the physical parameters related to the
76//model parameters. See the reference at the top of the
77//file for details
78//
79// this macro is currently only appicable to the monodisperse case and must be appropriately
80// modified before use
81//
82//Macro NumberDensity_FractalPolyCore()
83//     
84//      Variable nden,phi,r0,Df,corr,s0,vpoly,i0,rg
85//     
86//      if(Exists("coef_fraPolyCore")!=1)
87//              abort "You need to plot the unsmeared model first to create the coefficient table"
88//      Endif
89//     
90//      phi = coef_fraPolyCore[0]   // volume fraction of building blocks
91//      r0 = coef_fraPolyCore[1]    // building block radius
92//      Df = coef_fraPolyCore[2]    // fractal dimension
93//      corr = coef_fraPolyCore[3]  // fractal correlation length (of cluster)
94//     
95//      Print "mean building block radius (A) = ",r0
96//      Print "volume fraction = ",phi
97//     
98// //  average particle volume         
99//      vpoly = 4*Pi/3*r0^3
100//      nden = phi/vpoly                //nden in 1/A^3
101//      i0 = 1.0e8*phi*vpoly*(coef_fraPolyCore[4]-coef_fraPolyCore[5])^2  // 1/cm/sr
102//      rg = corr*( Df*(Df+1)/2 )^0.5
103//      s0 = exp(gammln(Df+1))*(corr/r0)^Df
104//      Print "number density (A^-3) = ",nden
105//      Print "Guinier radius (A) = ",rg
106//      Print "Aggregation number G = ",s0
107//      Print "Forward cross section of building blocks (cm-1 sr-1) I(0) = ",i0
108//      Print "Forward cross section of clusters (cm-1 sr-1) I(0) = ",i0*s0
109//End
110
111//AAO version, uses XOP if available
112// simply calls the original single point calculation with
113// a wave assignment (this will behave nicely if given point ranges)
114Function FractalPolyCore(cw,yw,xw) : FitFunc
115        Wave cw,yw,xw
116       
117#if exists("FractalPolyCoreX")
118        yw = FractalPolyCoreX(cw,xw)
119#else
120        yw = fFractalPolyCore(cw,xw)
121#endif
122        return(0)
123End
124
125//fractal scattering function
126Function fFractalPolyCore(w,x) :FitFunc
127        wave w
128        variable x
129       
130        variable r0,Df,corr,phi,sldp,sldm,bkg
131        variable pq,sq,ans,pd,thick,slds
132        phi=w[0]   // volume fraction of building block spheres...
133        r0=w[1]    //  radius of building block
134        pd = w[2]               // polydispersity of core
135        thick = w[3] // thell thickness
136        Df=w[4]     //  fractal dimension
137        corr=w[5] //  correlation length of fractal-like aggregates
138        sldp = w[6] // SLD of building block
139        slds = w[7]     //SLD of shell
140        sldm = w[8] // SLD of matrix or solution
141        bkg=w[9]  //  flat background
142       
143        //calculate P(q) for the spherical subunits, units cm-1 sr-1
144        Make/O/D/N=8 tmp_PolyCor
145        tmp_PolyCor[0] = phi
146        tmp_PolyCor[1] = r0
147        tmp_PolyCor[2] = pd
148        tmp_PolyCor[3] = thick
149        tmp_PolyCor[4] = sldp
150        tmp_PolyCor[5] = slds
151        tmp_PolyCor[6] = sldm
152        tmp_PolyCor[7] = 0
153
154#if exists("PolyCoreFormX")
155        pq = PolyCoreFormX(tmp_PolyCor,x)
156#else
157        pq = fPolyCoreForm(tmp_PolyCor,x)
158#endif
159
160        //calculate S(q)
161        sq = Df*exp(gammln(Df-1))*sin((Df-1)*atan(x*corr))
162        sq /= (x*r0)^Df * (1 + 1/(x*corr)^2)^((Df-1)/2)
163        sq += 1
164        //combine and return
165        ans = pq*sq + bkg
166
167        return (ans)
168End
169
170//wrapper to calculate the smeared model as an AAO-Struct
171// fills the struct and calls the ususal function with the STRUCT parameter
172//
173// used only for the dependency, not for fitting
174//
175Function fSmearedFractalPolyCore(coefW,yW,xW)
176        Wave coefW,yW,xW
177       
178        String str = getWavesDataFolder(yW,0)
179        String DF="root:"+str+":"
180       
181        WAVE resW = $(DF+str+"_res")
182       
183        STRUCT ResSmearAAOStruct fs
184        WAVE fs.coefW = coefW   
185        WAVE fs.yW = yW
186        WAVE fs.xW = xW
187        WAVE fs.resW = resW
188       
189        Variable err
190        err = SmearedFractalPolyCore(fs)
191       
192        return (0)
193End
194
195//the smeared model calculation
196Function SmearedFractalPolyCore(s) :FitFunc
197        Struct ResSmearAAOStruct &s
198
199//      the name of your unsmeared model (AAO) is the first argument
200        Smear_Model_20(FractalPolyCore,s.coefW,s.xW,s.yW,s.resW)
201
202        return(0)
203End
Note: See TracBrowser for help on using the repository browser.