source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/BimodalSchulzSpheres.ipf @ 151

Last change on this file since 151 was 151, checked in by srkline, 15 years ago

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

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