source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/Beaucage.ipf @ 61

Last change on this file since 61 was 61, checked in by srkline, 16 years ago

added a scale factor to the Beaucage model so that it can be properly used with Global Fitting (so that either USANS or SANS absolute scale can be tweaked as needed)
File size: 11.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3/////////////////////////////////////////////////////
4//
5// Plot's Greg Beaucage's Rg-power Law "model" of scattering
6// somewhat useful for identifying length scales, but short on
7// physical inerpretation of the real structure of the sample.
8//
9// up to 4 "levels" can be calculated
10// best to start with single level, and fit a small range of
11// the data, and add more levels as needed
12//
13// see the help file for the original references
14//
15Proc PlotBeau_One(num,qmin,qmax)
16        Variable num=256,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_b1,ywave_b1
22        xwave_b1 = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        make/o/d coef_b1 = {1,3,21,6e-4,2,0}
24        make/o/t parameters_b1 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","bkg (cm-1 sr-1)"}
25        Edit parameters_b1,coef_b1
26        ywave_b1 := OneLevel(coef_b1,xwave_b1)
27        Display ywave_b1 vs xwave_b1
28        ModifyGraph log=1,marker=29,msize=2,mode=4
29        Label bottom "q (\\S-1\\M)"
30        Label left "Intensity (cm\\S-1\\M)"
31        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
32End
33
34Proc PlotBeau_Two(num,qmin,qmax)
35        Variable num=256,qmin=0.001,qmax=0.7
36        Prompt num "Enter number of data points for model: "
37        Prompt qmin "Enter minimum q-value (^-1) for model: "
38        Prompt qmax "Enter maximum q-value (^-1) for model: "
39       
40        make/o/d/n=(num) xwave_b2,ywave_b2
41        xwave_b2 = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
42        make/o/d coef_b2 = {1,400,200,5e-6,4,3,21,6e-4,2,0}
43        make/o/t parameters_b2 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","bkg (cm-1 sr-1)"}
44        Edit parameters_b2,coef_b2
45        ywave_b2 := TwoLevel(coef_b2,xwave_b2)
46        Display ywave_b2 vs xwave_b2
47        ModifyGraph log=1,marker=29,msize=2,mode=4
48        Label bottom "q (\\S-1\\M)"
49        Label left "Intensity (cm\\S-1\\M)"
50        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
51End
52
53Proc PlotBeau_Three(num,qmin,qmax)
54        Variable num=256,qmin=0.001,qmax=0.7
55        Prompt num "Enter number of data points for model: "
56        Prompt qmin "Enter minimum q-value (^-1) for model: "
57        Prompt qmax "Enter maximum q-value (^-1) for model: "
58       
59        make/o/d/n=(num) xwave_b3,ywave_b3
60        xwave_b3 = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
61        make/o/d coef_b3 = {1,4000,600,2e-7,4,400,200,5e-6,4,3,21,6e-4,2,0}
62        make/o/t parameters_b3 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","G3 (cm-1 sr-1)","Rg3  (A)","B3 (cm-1 sr-1)","Pow3","bkg (cm-1)"}
63        Edit parameters_b3,coef_b3
64        ywave_b3 := ThreeLevel(coef_b3,xwave_b3)
65        Display ywave_b3 vs xwave_b3
66        ModifyGraph log=1,marker=29,msize=2,mode=4
67        Label bottom "q (\\S-1\\M)"
68        Label left "Intensity (cm\\S-1\\M)"
69        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
70End
71
72Proc PlotBeau_Four(num,qmin,qmax)
73        Variable num=256,qmin=0.001,qmax=0.7
74        Prompt num "Enter number of data points for model: "
75        Prompt qmin "Enter minimum q-value (^-1) for model: "
76        Prompt qmax "Enter maximum q-value (^-1) for model: "
77       
78        make/o/d/n=(num) xwave_b4,ywave_b4
79        xwave_b4 = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
80        make/o/d coef_b4 = {1,40000,2000,1e-8,4,4000,600,2e-7,4,400,200,5e-6,4,3,21,6e-4,2,0}
81        make/o/t parameters_b4 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","G3 (cm-1 sr-1)","Rg3  (A)","B3 (cm-1 sr-1)","Pow3","G4 (cm-1 sr-1)","Rg4  (A)","B4 (cm-1 sr-1)","Pow4","bkg (cm-1)"}
82        Edit parameters_b4,coef_b4
83        ywave_b4 := FourLevel(coef_b4,xwave_b4)
84        Display ywave_b4 vs xwave_b4
85        ModifyGraph log=1,marker=29,msize=2,mode=4
86        Label bottom "q (\\S-1\\M)"
87        Label left "Intensity (cm\\S-1\\M)"
88        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
89End
90
91/////////// macros for smeared model calculations
92
93Proc PlotSmearedBeau_OneLevel()                                                         
94        //no input parameters necessary, it MUST use the experimental q-values
95        // from the experimental data read in from an AVE/QSIG data file
96        If(ResolutionWavesMissing())            //part of GaussUtils
97                Abort
98        endif
99       
100        // Setup parameter table for model function
101        Make/O/D smear_coef_b1 ={1,3,21,6e-4,2,0}                                       
102        make/o/t smear_parameters_b1 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","bkg (cm-1 sr-1)"} 
103        Edit smear_parameters_b1,smear_coef_b1                                 
104       
105        // output smeared intensity wave, dimensions are identical to experimental QSIG values
106        // make extra copy of experimental q-values for easy plotting
107        Duplicate/O $gQvals smeared_b1,smeared_qvals                           
108        SetScale d,0,0,"1/cm",smeared_b1                                                       
109
110        smeared_b1 := SmearedOneLevel(smear_coef_b1,$gQvals)           
111        Display smeared_b1 vs smeared_qvals                                                                     
112        ModifyGraph log=1,marker=29,msize=2,mode=4
113        Label bottom "q (\\S-1\\M)"
114        Label left "Intensity (cm\\S-1\\M)"
115        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
116End
117
118Proc PlotSmearedBeau_TwoLevel()                                                         
119        //no input parameters necessary, it MUST use the experimental q-values
120        // from the experimental data read in from an AVE/QSIG data file
121        If(ResolutionWavesMissing())            //part of GaussUtils
122                Abort
123        endif
124       
125        // Setup parameter table for model function
126        Make/O/D smear_coef_b2 = {1,400,200,5e-6,4,3,21,6e-4,2,0}                               
127        make/o/t smear_parameters_b2 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","bkg (cm-1 sr-1)"}     
128        Edit smear_parameters_b2,smear_coef_b2                                 
129       
130        // output smeared intensity wave, dimensions are identical to experimental QSIG values
131        // make extra copy of experimental q-values for easy plotting
132        Duplicate/O $gQvals smeared_b2,smeared_qvals                           
133        SetScale d,0,0,"1/cm",smeared_b2                                                       
134
135        smeared_b2 := SmearedTwoLevel(smear_coef_b2,$gQvals)           
136        Display smeared_b2 vs smeared_qvals                                                                     
137        ModifyGraph log=1,marker=29,msize=2,mode=4
138        Label bottom "q (\\S-1\\M)"
139        Label left "Intensity (cm\\S-1\\M)"
140        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
141End
142
143Proc PlotSmearedBeau_ThreeLevel()                                                               
144        //no input parameters necessary, it MUST use the experimental q-values
145        // from the experimental data read in from an AVE/QSIG data file
146        If(ResolutionWavesMissing())            //part of GaussUtils
147                Abort
148        endif
149       
150        // Setup parameter table for model function
151        Make/O/D smear_coef_b3 = {1,4000,600,2e-7,4,400,200,5e-6,4,3,21,6e-4,2,0}
152        make/o/t smear_parameters_b3 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","G3 (cm-1 sr-1)","Rg3  (A)","B3 (cm-1 sr-1)","Pow3","bkg (cm-1)"}
153        Edit smear_parameters_b3,smear_coef_b3                                 
154       
155        // output smeared intensity wave, dimensions are identical to experimental QSIG values
156        // make extra copy of experimental q-values for easy plotting
157        Duplicate/O $gQvals smeared_b3,smeared_qvals                           
158        SetScale d,0,0,"1/cm",smeared_b3                                                       
159
160        smeared_b3 := SmearedThreeLevel(smear_coef_b3,$gQvals)         
161        Display smeared_b3 vs smeared_qvals                                                                     
162        ModifyGraph log=1,marker=29,msize=2,mode=4
163        Label bottom "q (\\S-1\\M)"
164        Label left "Intensity (cm\\S-1\\M)"
165        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
166End
167
168Proc PlotSmearedBeau_FourLevel()                                                               
169        //no input parameters necessary, it MUST use the experimental q-values
170        // from the experimental data read in from an AVE/QSIG data file
171        If(ResolutionWavesMissing())            //part of GaussUtils
172                Abort
173        endif
174       
175        // Setup parameter table for model function
176        Make/O/D smear_coef_b4 = {1,40000,2000,1e-8,4,4000,600,2e-7,4,400,200,5e-6,4,3,21,6e-4,2,0}
177        Make/o/t smear_parameters_b4 = {"scale","G1 (cm-1 sr-1)","Rg1  (A)","B1 (cm-1 sr-1)","Pow1","G2 (cm-1 sr-1)","Rg2  (A)","B2 (cm-1 sr-1)","Pow2","G3 (cm-1 sr-1)","Rg3  (A)","B3 (cm-1 sr-1)","Pow3","G4 (cm-1 sr-1)","Rg4  (A)","B4 (cm-1 sr-1)","Pow4","bkg (cm-1)"}
178        Edit smear_parameters_b4,smear_coef_b4                                 
179       
180        // output smeared intensity wave, dimensions are identical to experimental QSIG values
181        // make extra copy of experimental q-values for easy plotting
182        Duplicate/O $gQvals smeared_b4,smeared_qvals                           
183        SetScale d,0,0,"1/cm",smeared_b4                                                       
184
185        smeared_b4 := SmearedFourLevel(smear_coef_b4,$gQvals)           
186        Display smeared_b4 vs smeared_qvals                                                                     
187        ModifyGraph log=1,marker=29,msize=2,mode=4
188        Label bottom "q (\\S-1\\M)"
189        Label left "Intensity (cm\\S-1\\M)"
190        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
191End
192
193//////////Function definitions
194
195Function OneLevel(w,x) :FitFunc
196        Wave w
197        Variable x
198       
199        Variable ans,erf1,prec=1e-15
200        Variable G1,Rg1,B1,Pow1,bkg
201       
202        G1 = w[0]
203        Rg1 = w[1]
204        B1 = w[2]
205        Pow1 = w[3]
206        bkg = w[4]
207       
208        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
209       
210        ans = G1*exp(-x*x*Rg1*Rg1/3)
211        ans += B1*(erf1^3/x)^Pow1
212        ans += bkg
213       
214        Return (ans)
215End
216
217Function TwoLevel(w,x) :FitFunc
218        Wave w
219        Variable x
220       
221        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
222        Variable erf1,erf2,prec=1e-15
223       
224        //Rsub = Rs
225        G1 = w[0]       //equivalent to I(0)
226        Rg1 = w[1]
227        B1 = w[2]
228        Pow1 = w[3]
229        G2 = w[4]
230        Rg2 = w[5]
231        B2 = w[6]
232        Pow2 = w[7]
233        bkg = w[8]
234       
235        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
236        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
237        //Print erf1
238       
239        ans = G1*exp(-x*x*Rg1*Rg1/3)
240        ans += B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
241        ans += G2*exp(-x*x*Rg2*Rg2/3)
242        ans += B2*(erf2^3/x)^Pow2
243       
244        ans += bkg
245       
246        Return(ans)
247End
248
249
250Function ThreeLevel(w,x) :FitFunc
251        Wave w
252        Variable x
253       
254        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
255        Variable G3,Rg3,B3,Pow3,erf3
256        Variable erf1,erf2,prec=1e-15
257       
258        //Rsub = Rs
259        G1 = w[0]       //equivalent to I(0)
260        Rg1 = w[1]
261        B1 = w[2]
262        Pow1 = w[3]
263        G2 = w[4]
264        Rg2 = w[5]
265        B2 = w[6]
266        Pow2 = w[7]
267        G3 = w[8]
268        Rg3 = w[9]
269        B3 = w[10]
270        Pow3 = w[11]
271        bkg = w[12]
272       
273        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
274        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
275        erf3 = erf( (x*Rg3/sqrt(6)) ,prec)
276        //Print erf1
277       
278        ans = G1*exp(-x*x*Rg1*Rg1/3) + B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
279        ans += G2*exp(-x*x*Rg2*Rg2/3) + B2*exp(-x*x*Rg3*Rg3/3)*(erf2^3/x)^Pow2
280        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*(erf3^3/x)^Pow3
281       
282        ans += bkg
283       
284        Return(ans)
285End
286
287Function FourLevel(w,x) :FitFunc
288        Wave w
289        Variable x
290       
291        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
292        Variable G3,Rg3,B3,Pow3,erf3
293        Variable G4,Rg4,B4,Pow4,erf4
294        Variable erf1,erf2,prec=1e-15
295       
296        //Rsub = Rs
297        G1 = w[0]       //equivalent to I(0)
298        Rg1 = w[1]
299        B1 = w[2]
300        Pow1 = w[3]
301        G2 = w[4]
302        Rg2 = w[5]
303        B2 = w[6]
304        Pow2 = w[7]
305        G3 = w[8]
306        Rg3 = w[9]
307        B3 = w[10]
308        Pow3 = w[11]
309        G4 = w[12]
310        Rg4 = w[13]
311        B4 = w[14]
312        Pow4 = w[15]
313        bkg = w[16]
314       
315        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
316        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
317        erf3 = erf( (x*Rg3/sqrt(6)) ,prec)
318        erf4 = erf( (x*Rg4/sqrt(6)) ,prec)
319       
320        ans = G1*exp(-x*x*Rg1*Rg1/3) + B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
321        ans += G2*exp(-x*x*Rg2*Rg2/3) + B2*exp(-x*x*Rg3*Rg3/3)*(erf2^3/x)^Pow2
322        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*exp(-x*x*Rg4*Rg4/3)*(erf3^3/x)^Pow3
323        ans += G4*exp(-x*x*Rg4*Rg4/3) + B4*(erf4^3/x)^Pow4
324       
325        ans += bkg
326       
327        Return(ans)
328End
329
330Function SmearedOneLevel(w,x) :FitFunc
331        Wave w
332        Variable x
333       
334//      Variable timer=StartMSTimer
335        Variable ans
336        SVAR sq = gSig_Q
337        SVAR qb = gQ_bar
338        SVAR sh = gShadow
339        SVAR gQ = gQVals
340       
341        ans = Smear_Model_20(OneLevel,$sq,$qb,$sh,$gQ,w,x)     
342
343//      Print "HS elapsed time(s) = ",StopMSTimer(timer)*1e-6
344        return(ans)
345End
346
347Function SmearedTwoLevel(w,x) :FitFunc
348        Wave w
349        Variable x
350       
351        Variable ans
352        SVAR sq = gSig_Q
353        SVAR qb = gQ_bar
354        SVAR sh = gShadow
355        SVAR gQ = gQVals
356       
357        ans = Smear_Model_20(TwoLevel,$sq,$qb,$sh,$gQ,w,x)     
358
359        return(ans)
360End
361
362Function SmearedThreeLevel(w,x) :FitFunc
363        Wave w
364        Variable x
365       
366        Variable ans
367        SVAR sq = gSig_Q
368        SVAR qb = gQ_bar
369        SVAR sh = gShadow
370        SVAR gQ = gQVals
371       
372        ans = Smear_Model_20(ThreeLevel,$sq,$qb,$sh,$gQ,w,x)   
373
374        return(ans)
375End
376
377Function SmearedFourLevel(w,x) :FitFunc
378        Wave w
379        Variable x
380       
381        Variable ans
382        SVAR sq = gSig_Q
383        SVAR qb = gQ_bar
384        SVAR sh = gShadow
385        SVAR gQ = gQVals
386       
387        ans = Smear_Model_20(FourLevel,$sq,$qb,$sh,$gQ,w,x)     
388
389        return(ans)
390End
Note: See TracBrowser for help on using the repository browser.