source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/NewModels_2006/BimodalSchulzSpheres.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 7.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// this model calculation is for the scattered intensity from a dispersion of polydisperse spheres
5// hard sphere interactions are NOT included
6// the polydispersity in radius is a Schulz distribution
7//
8// TWO polulations of spheres are considered
9//
10// 31 DEC 03 SRK
11////////////////////////////////////////////////
12#include "SchulzSpheres"
13
14Proc PlotBimodalSchulzSpheres(num,qmin,qmax)
15        Variable num=128,qmin=0.001,qmax=0.7
16        Prompt num "Enter number of data points for model: "
17        Prompt qmin "Enter minimum q-value (A^-1) for model: "
18        Prompt qmax "Enter maximum q-value (A^-1) for model: "
19       
20        Make/O/D/n=(num) xwave_bss,ywave_bss
21        xwave_bss =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
22        Make/O/D coef_bss = {0.01,200,0.2,1e-6,0.05,25,0.2,1e-6,6.4e-6,0.001}
23        make/o/t/N=10 parameters_bss
24        parameters_bss[0,3] = {"volume fraction(1)","Radius (1) (A)","polydispersity(1)","SLD(1) (A^-2)"}
25        parameters_bss[4,9] = {"volume fraction(2)","Radius (2)","polydispersity(2)","SLD(2)","SLD (solvent)","background (cm-1 sr-1)"}
26        Edit parameters_bss,coef_bss
27        ywave_bss := BimodalSchulzSpheres(coef_bss,xwave_bss)
28        Display ywave_bss vs xwave_bss
29        ModifyGraph log=1,marker=29,msize=2,mode=4
30        Label bottom "q (A\\S-1\\M)"
31        Label left "Intensity (cm\\S-1\\M)"
32        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
33End
34
35
36///////////////////////////////////////////////////////////
37
38Proc PlotSmearedBimodalSchulzSpheres()         
39        //no input parameters necessary, it MUST use the experimental q-values
40        // from the experimental data read in from an AVE/QSIG data file
41       
42        // if no gQvals wave, data must not have been loaded => abort
43        if(ResolutionWavesMissing())
44                Abort
45        endif
46       
47        // Setup parameter table for model function
48        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}
49        make/o/t/N=10 smear_parameters_bss
50        smear_parameters_bss[0,3] = {"volume fraction(1)","Radius (1) (A)","polydispersity(1)","SLD(1) (A^-2)"}
51        smear_parameters_bss[4,9] = {"volume fraction(2)","Radius (2)","polydispersity(2)","SLD(2)","SLD (solvent)","background (cm-1 sr-1)"}
52        Edit smear_parameters_bss,smear_coef_bss
53       
54        // output smeared intensity wave, dimensions are identical to experimental QSIG values
55        // make extra copy of experimental q-values for easy plotting
56        Duplicate/O $gQvals smeared_bss,smeared_qvals
57        SetScale d,0,0,"1/cm",smeared_bss
58
59        smeared_bss := SmearedBimodalSchulzSpheres(smear_coef_bss,$gQvals)
60        Display smeared_bss vs smeared_qvals
61        ModifyGraph log=1,marker=29,msize=2,mode=4
62        Label bottom "q (A\\S-1\\M)"
63        Label left "Intensity (cm\\S-1\\M)"
64        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
65End
66
67//  Calculates some characteristic parameters for bimodal Shulz distribution
68Macro NumberDensity_Bimodal()
69
70        Variable nden1,nden2,phi1,phi2,R1,R2,Ravg,p1,p2,Rg1,Rg2,I1_0,I2_0,I0,Sv1,Sv2,Sv,vpoly1,vpoly2
71        Variable z1,z2,v2poly1,v2poly2
72        if(WaveExists(coef_bss)==0)
73                abort "You need to plot the model first to create the coefficient table"
74        Endif
75       
76        phi1 = coef_bss[0]  // volume fraction, mode 1
77        phi2 = coef_bss[4]  // volume fraction, mode 1
78        R1 = coef_bss[1]  // mean radius, mode 1(A)
79        R2 = coef_bss[5]  // mean radius, mode 1(A)
80        p1 = coef_bss[2]  // polydispersity, mode 1
81        p2 = coef_bss[6]  // polydispersity, mode 1
82                       
83        z1 = (1/p1)^2-1
84        z2 = (1/p2)^2-1
85        //  average particle volume     
86        vpoly1 = 4*Pi/3*(z1+3)*(z1+2)/(z1+1)/(z1+1)*r1^3
87        vpoly2 = 4*Pi/3*(z2+3)*(z2+2)/(z2+1)/(z2+1)*r2^3
88        //average particle volume^2     
89        v2poly1 = (4*Pi/3)^2*(z1+6)*(z1+5)*(z1+4)*(z1+3)*(z1+2)/((z1+1)^5)*r1^6
90        v2poly2 = (4*Pi/3)^2*(z2+6)*(z2+5)*(z2+4)*(z2+3)*(z2+2)/((z2+1)^5)*r2^6
91        nden1 = phi1/vpoly1             //nden in 1/A^3
92        nden2 = phi2/vpoly2             //nden in 1/A^3
93
94        rg1 = r1*((3*(z1+8)*(z1+7))/5/(z1+1)/(z1+1))^0.5   // in A
95        rg2 = r2*((3*(z2+8)*(z2+7))/5/(z2+1)/(z2+1))^0.5   // in A
96        sv1 = 1.0e8*3*phi1*(z1+1)/R1/(z1+3)  // in 1/cm
97        sv2 = 1.0e8*3*phi2*(z2+1)/R2/(z2+3)  // in 1/cm
98        I1_0 = 1.0e8*nden1*v2poly1*(coef_bss[3]-coef_bss[8])^2  // 1/cm/sr
99        I2_0 = 1.0e8*nden2*v2poly2*(coef_bss[7]-coef_bss[8])^2  // 1/cm/sr
100       
101        Print "mode 1 number density (A^-3) = ",nden1
102        Print "mode 2 number density (A^-3) = ",nden2
103
104        Ravg = (nden1*R1+nden2*R2)/(nden1+nden2)
105
106        Print "mean radius, mode 1 (A) = ",R1
107        Print "mean radius, mode 2 (A) = ",R2
108        Print "mean radius, total     (A) = ",Ravg
109        Print "polydispersity, mode 1 (sig/avg) = ",p1
110        Print "polydispersity, mode 2 (sig/avg) = ",p2
111        Print "volume fraction, mode 1 = ",phi1
112        Print "volume fraction, mode 2 = ",phi2
113
114        Print "Guinier Radius, mode 1 (A) = ",Rg1
115        Print "Guinier Radius, mode 2 (A) = ",Rg2
116        I0 = I1_0+I2_0
117        Print "Forward scattering cross-section, mode 1 (cm-1 sr-1) I(0)= ",I1_0
118        Print "Forward scattering cross-section, mode 2 (cm-1 sr-1) I(0)= ",I2_0
119        Print "Forward scattering cross-section, total     (cm-1 sr-1) I(0)= ",I0
120        Sv = Sv1+Sv2
121        Print "Interfacial surface area per unit sample volume, mode 1 (cm-1) Sv= ",Sv1
122        Print "Interfacial surface area per unit sample volume, mode 2 (cm-1) Sv= ",Sv2
123        Print "Interfacial surface area per unit sample volume, total     (cm-1) Sv= ",Sv       
124End
125
126
127//   Plots bimodal size distribution
128Macro Plot_Bimodal_Distribution()
129
130        variable p1,p2,r1,r2,z1,z2,phi1,phi2,f1,f2,nden1,nden2,vpoly1,vpoly2,maxr
131       
132        if(WaveExists(coef_bss)==0)
133                abort "You need to plot the model first to create the coefficient table"
134        Endif
135        phi1 = coef_bss[0]  // volume fraction, mode 1
136        phi2 = coef_bss[4]  // volume fraction, mode 1
137        R1 = coef_bss[1]  // mean radius, mode 1(A)
138        R2 = coef_bss[5]  // mean radius, mode 1(A)
139        p1 = coef_bss[2]  // polydispersity, mode 1
140        p2 = coef_bss[6]  // polydispersity, mode 1
141                       
142        z1 = (1/p1)^2-1
143        z2 = (1/p2)^2-1
144        //  average particle volume     
145        vpoly1 = 4*Pi/3*(z1+3)*(z1+2)/(z1+1)/(z1+1)*r1^3
146        vpoly2 = 4*Pi/3*(z2+3)*(z2+2)/(z2+1)/(z2+1)*r2^3
147
148        nden1 = phi1/vpoly1             //nden in 1/A^3
149        nden2 = phi2/vpoly2             //nden in 1/A^3
150        f1 = nden1/(nden1+nden2)
151        f2 = nden2/(nden1+nden2)
152
153        Make/O/D/N=1000 Bimodal_Schulz_distribution
154        if (r1>r2) then
155           maxr =  r1*(1+6*p1)
156        else
157           maxr =  r2*(1+6*p2)
158        endif
159
160        SetScale/I x, 0, maxr, Bimodal_Schulz_distribution
161        Bimodal_Schulz_distribution = f1*Schulz_Point_bss(x,r1,z1)+f2*Schulz_Point_bss(x,r2,z2)
162        Display Bimodal_Schulz_distribution
163        Label left "f(R) (normalized)"
164        Label bottom "R (A)"
165        legend
166End
167
168///////////////////////////////////////////////////////////////
169// unsmeared model calculation
170///////////////////////////
171Function BimodalSchulzSpheres(w,k) : FitFunc
172        Wave w                  // the coefficient wave
173        Variable k              // the x values, as a variable
174
175        Variable ans=0
176        Make/O/D/N=6 temp_coef_1,temp_coef_2            //coefficient waves for each population
177        temp_coef_1[0] = w[0]
178        temp_coef_1[1] = w[1]
179        temp_coef_1[2] = w[2]
180        temp_coef_1[3] = w[3]
181        temp_coef_1[4] = w[8]
182        temp_coef_1[5] = 0
183
184//second population
185        temp_coef_2[0] = w[4]
186        temp_coef_2[1] = w[5]
187        temp_coef_2[2] = w[6]
188        temp_coef_2[3] = w[7]
189        temp_coef_2[4] = w[8]
190        temp_coef_2[5] = 0              //always zero - background is added in the final step
191       
192        //calculate both models and sum (add background here)
193        ans = SchulzSpheres(temp_coef_1,k)
194        ans += SchulzSpheres(temp_coef_2,k)
195        ans += w[9]             //background
196       
197        return(ans)
198End
199
200Function Schulz_Point_bss(x,avg,zz)
201        Variable x,avg,zz
202       
203        Variable dr
204       
205        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
206        return (exp(dr))
207End
208
209// this is all there is to the smeared calculation!
210Function SmearedBimodalSchulzSpheres(w,x) :FitFunc
211        Wave w
212        Variable x
213       
214        Variable ans
215        SVAR sq = gSig_Q
216        SVAR qb = gQ_bar
217        SVAR sh = gShadow
218        SVAR gQ = gQVals
219       
220        //the name of your unsmeared model is the first argument
221        ans = Smear_Model_20(BimodalSchulzSpheres,$sq,$qb,$sh,$gQ,w,x)
222
223        return(ans)
224End
225
Note: See TracBrowser for help on using the repository browser.