source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/BimodalSchulzSpheres_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

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       
37        AddModelToStrings("BimodalSchulzSpheres","coef_bss","parameters_bss","bss")
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:
76        AddModelToStrings("SmearedBimodalSchulzSpheres","smear_coef_bss","smear_parameters_bss","bss")
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
Note: See TracBrowser for help on using the repository browser.