# source:sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/BimodalSchulzSpheres_v40.ipf@379

Last change on this file since 379 was 379, checked in by srkline, 14 years ago

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

File size: 8.5 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// this model calculation is for the scattered intensity from a dispersion of polydisperse spheres
6// hard sphere interactions are NOT included
7// the polydispersity in radius is a Schulz distribution
8//
9// TWO polulations of spheres are considered
10//
11// 31 DEC 03 SRK
12////////////////////////////////////////////////
13#include "SchulzSpheres_v40"
14
15Proc PlotBimodalSchulzSpheres(num,qmin,qmax)
16        Variable num=128,qmin=0.001,qmax=0.7
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_bss,ywave_bss
22        xwave_bss =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        Make/O/D coef_bss = {0.01,200,0.2,1e-6,0.05,25,0.2,1e-6,6.4e-6,0.001}
24        make/o/t/N=10 parameters_bss
25        parameters_bss[0,3] = {"volume fraction(1)","Radius (1) (A)","polydispersity(1)","SLD(1) (A^-2)"}
26        parameters_bss[4,9] = {"volume fraction(2)","Radius (2)","polydispersity(2)","SLD(2)","SLD (solvent)","background (cm-1 sr-1)"}
27        Edit parameters_bss,coef_bss
28
29        Variable/G root:g_bss
30        g_bss := BimodalSchulzSpheres(coef_bss,ywave_bss,xwave_bss)
31        Display ywave_bss vs xwave_bss
32        ModifyGraph log=1,marker=29,msize=2,mode=4
33        Label bottom "q (A\\S-1\\M)"
34        Label left "Intensity (cm\\S-1\\M)"
35        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
36
38End
39
40
41///////////////////////////////////////////////////////////
42// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
43Proc PlotSmearedBimodalSchulzSpheres(str)
44        String str
45        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
46
47        // if any of the resolution waves are missing => abort
48        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
49                Abort
50        endif
51
52        SetDataFolder \$("root:"+str)
53
54        // Setup parameter table for model function
55        Make/O/D smear_coef_bss = {0.01,200,0.2,1e-6,0.05,25,0.2,1e-6,6.4e-6,0.001}
56        make/o/t/N=10 smear_parameters_bss
57        smear_parameters_bss[0,3] = {"volume fraction(1)","Radius (1) (A)","polydispersity(1)","SLD(1) (A^-2)"}
58        smear_parameters_bss[4,9] = {"volume fraction(2)","Radius (2)","polydispersity(2)","SLD(2)","SLD (solvent)","background (cm-1 sr-1)"}
59        Edit smear_parameters_bss,smear_coef_bss
60
61        // output smeared intensity wave, dimensions are identical to experimental QSIG values
62        // make extra copy of experimental q-values for easy plotting
63        Duplicate/O \$(str+"_q") smeared_bss,smeared_qvals
64        SetScale d,0,0,"1/cm",smeared_bss
65
66        Variable/G gs_bss=0
67        gs_bss := fSmearedBimodalSchulzSpheres(smear_coef_bss,smeared_bss,smeared_qvals)        //this wrapper fills the STRUCT
68
69        Display smeared_bss vs smeared_qvals
70        ModifyGraph log=1,marker=29,msize=2,mode=4
71        Label bottom "q (A\\S-1\\M)"
72        Label left "Intensity (cm\\S-1\\M)"
73        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
74
75        SetDataFolder root:
77End
78
79
80//  Calculates some characteristic parameters for bimodal Shulz distribution
81Macro NumberDensity_Bimodal()
82
83        Variable nden1,nden2,phi1,phi2,R1,R2,Ravg,p1,p2,Rg1,Rg2,I1_0,I2_0,I0,Sv1,Sv2,Sv,vpoly1,vpoly2
84        Variable z1,z2,v2poly1,v2poly2
85
86        if(Exists("coef_bss")!=1)
87                abort "You need to plot the unsmeared model first to create the coefficient table"
88        Endif
89
90        phi1 = coef_bss[0]  // volume fraction, mode 1
91        phi2 = coef_bss[4]  // volume fraction, mode 1
92        R1 = coef_bss[1]  // mean radius, mode 1(A)
93        R2 = coef_bss[5]  // mean radius, mode 1(A)
94        p1 = coef_bss[2]  // polydispersity, mode 1
95        p2 = coef_bss[6]  // polydispersity, mode 1
96
97        z1 = (1/p1)^2-1
98        z2 = (1/p2)^2-1
99        //  average particle volume
100        vpoly1 = 4*Pi/3*(z1+3)*(z1+2)/(z1+1)/(z1+1)*r1^3
101        vpoly2 = 4*Pi/3*(z2+3)*(z2+2)/(z2+1)/(z2+1)*r2^3
102        //average particle volume^2
103        v2poly1 = (4*Pi/3)^2*(z1+6)*(z1+5)*(z1+4)*(z1+3)*(z1+2)/((z1+1)^5)*r1^6
104        v2poly2 = (4*Pi/3)^2*(z2+6)*(z2+5)*(z2+4)*(z2+3)*(z2+2)/((z2+1)^5)*r2^6
105        nden1 = phi1/vpoly1             //nden in 1/A^3
106        nden2 = phi2/vpoly2             //nden in 1/A^3
107
108        rg1 = r1*((3*(z1+8)*(z1+7))/5/(z1+1)/(z1+1))^0.5   // in A
109        rg2 = r2*((3*(z2+8)*(z2+7))/5/(z2+1)/(z2+1))^0.5   // in A
110        sv1 = 1.0e8*3*phi1*(z1+1)/R1/(z1+3)  // in 1/cm
111        sv2 = 1.0e8*3*phi2*(z2+1)/R2/(z2+3)  // in 1/cm
112        I1_0 = 1.0e8*nden1*v2poly1*(coef_bss[3]-coef_bss[8])^2  // 1/cm/sr
113        I2_0 = 1.0e8*nden2*v2poly2*(coef_bss[7]-coef_bss[8])^2  // 1/cm/sr
114
115        Print "mode 1 number density (A^-3) = ",nden1
116        Print "mode 2 number density (A^-3) = ",nden2
117
118        Ravg = (nden1*R1+nden2*R2)/(nden1+nden2)
119
120        Print "mean radius, mode 1 (A) = ",R1
121        Print "mean radius, mode 2 (A) = ",R2
122        Print "mean radius, total     (A) = ",Ravg
123        Print "polydispersity, mode 1 (sig/avg) = ",p1
124        Print "polydispersity, mode 2 (sig/avg) = ",p2
125        Print "volume fraction, mode 1 = ",phi1
126        Print "volume fraction, mode 2 = ",phi2
127
128        Print "Guinier Radius, mode 1 (A) = ",Rg1
129        Print "Guinier Radius, mode 2 (A) = ",Rg2
130        I0 = I1_0+I2_0
131        Print "Forward scattering cross-section, mode 1 (cm-1 sr-1) I(0)= ",I1_0
132        Print "Forward scattering cross-section, mode 2 (cm-1 sr-1) I(0)= ",I2_0
133        Print "Forward scattering cross-section, total     (cm-1 sr-1) I(0)= ",I0
134        Sv = Sv1+Sv2
135        Print "Interfacial surface area per unit sample volume, mode 1 (cm-1) Sv= ",Sv1
136        Print "Interfacial surface area per unit sample volume, mode 2 (cm-1) Sv= ",Sv2
137        Print "Interfacial surface area per unit sample volume, total     (cm-1) Sv= ",Sv
138End
139
140
141//   Plots bimodal size distribution
142Macro Plot_Bimodal_Distribution()
143
144        variable p1,p2,r1,r2,z1,z2,phi1,phi2,f1,f2,nden1,nden2,vpoly1,vpoly2,maxr
145
146        if(Exists("coef_bss")!=1)
147                abort "You need to plot the unsmeared model first to create the coefficient table"
148        Endif
149        phi1 = coef_bss[0]  // volume fraction, mode 1
150        phi2 = coef_bss[4]  // volume fraction, mode 1
151        R1 = coef_bss[1]  // mean radius, mode 1(A)
152        R2 = coef_bss[5]  // mean radius, mode 1(A)
153        p1 = coef_bss[2]  // polydispersity, mode 1
154        p2 = coef_bss[6]  // polydispersity, mode 1
155
156        z1 = (1/p1)^2-1
157        z2 = (1/p2)^2-1
158        //  average particle volume
159        vpoly1 = 4*Pi/3*(z1+3)*(z1+2)/(z1+1)/(z1+1)*r1^3
160        vpoly2 = 4*Pi/3*(z2+3)*(z2+2)/(z2+1)/(z2+1)*r2^3
161
162        nden1 = phi1/vpoly1             //nden in 1/A^3
163        nden2 = phi2/vpoly2             //nden in 1/A^3
164        f1 = nden1/(nden1+nden2)
165        f2 = nden2/(nden1+nden2)
166
167        Make/O/D/N=1000 Bimodal_Schulz_distribution
168        if (r1>r2) then
169           maxr =  r1*(1+6*p1)
170        else
171           maxr =  r2*(1+6*p2)
172        endif
173
174        SetScale/I x, 0, maxr, Bimodal_Schulz_distribution
175        Bimodal_Schulz_distribution = f1*Schulz_Point_bss(x,r1,z1)+f2*Schulz_Point_bss(x,r2,z2)
176        Display Bimodal_Schulz_distribution
177        Label left "f(R) (normalized)"
178        Label bottom "R (A)"
179        legend
180End
181
182///////////////////////////////////////////////////////////////
183// unsmeared model calculation
184///////////////////////////
185// now AAO function
186Function BimodalSchulzSpheres(w,yw,xw) : FitFunc
187        Wave w,yw,xw                    // the coefficient wave, y, x
188
189        Variable ans=0
190        Make/O/D/N=6 temp_coef_1,temp_coef_2            //coefficient waves for each population
191        temp_coef_1[0] = w[0]
192        temp_coef_1[1] = w[1]
193        temp_coef_1[2] = w[2]
194        temp_coef_1[3] = w[3]
195        temp_coef_1[4] = w[8]
196        temp_coef_1[5] = 0
197
198//second population
199        temp_coef_2[0] = w[4]
200        temp_coef_2[1] = w[5]
201        temp_coef_2[2] = w[6]
202        temp_coef_2[3] = w[7]
203        temp_coef_2[4] = w[8]
204        temp_coef_2[5] = 0              //always zero - background is added in the final step
205
206        //calculate both models and sum (add background here)
207        Duplicate/O xw tmp_ss1,tmp_ss2
208        SchulzSpheres(temp_coef_1,tmp_ss1,xw)
209        SchulzSpheres(temp_coef_2,tmp_ss2,xw)
210        yw = tmp_ss1 + tmp_ss2
211        yw += w[9]              //background
212
213        return(0)
214End
215
216Function Schulz_Point_bss(x,avg,zz)
217        Variable x,avg,zz
218
219        Variable dr
220
221        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
222        return (exp(dr))
223End
224
225//wrapper to calculate the smeared model as an AAO-Struct
226// fills the struct and calls the ususal function with the STRUCT parameter
227//
228// used only for the dependency, not for fitting
229//
230Function fSmearedBimodalSchulzSpheres(coefW,yW,xW)
231        Wave coefW,yW,xW
232
233        String str = getWavesDataFolder(yW,0)
234        String DF="root:"+str+":"
235
236        WAVE resW = \$(DF+str+"_res")
237
238        STRUCT ResSmearAAOStruct fs
239        WAVE fs.coefW = coefW
240        WAVE fs.yW = yW
241        WAVE fs.xW = xW
242        WAVE fs.resW = resW
243
244        Variable err
245        err = SmearedBimodalSchulzSpheres(fs)
246
247        return (0)
248End
249
250// this is all there is to the smeared calculation!
251Function SmearedBimodalSchulzSpheres(s) :FitFunc
252        Struct ResSmearAAOStruct &s
253
254//      the name of your unsmeared model (AAO) is the first argument
255        Smear_Model_20(BimodalSchulzSpheres,s.coefW,s.xW,s.yW,s.resW)
256
257        return(0)
258End
259
260
261
Note: See TracBrowser for help on using the repository browser.