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

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

scale factor added to final Beaucage model
File size: 11.7 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,scale
201       
202        scale = w[0]
203        G1 = w[1]
204        Rg1 = w[2]
205        B1 = w[3]
206        Pow1 = w[4]
207        bkg = w[5]
208       
209        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
210       
211        ans = G1*exp(-x*x*Rg1*Rg1/3)
212        ans += B1*(erf1^3/x)^Pow1
213        ans *= scale
214        ans += bkg
215       
216        Return (ans)
217End
218
219Function TwoLevel(w,x) :FitFunc
220        Wave w
221        Variable x
222       
223        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
224        Variable erf1,erf2,prec=1e-15,scale
225       
226        //Rsub = Rs
227        scale = w[0]
228        G1 = w[1]       //equivalent to I(0)
229        Rg1 = w[2]
230        B1 = w[3]
231        Pow1 = w[4]
232        G2 = w[5]
233        Rg2 = w[6]
234        B2 = w[7]
235        Pow2 = w[8]
236        bkg = w[9]
237       
238        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
239        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
240        //Print erf1
241       
242        ans = G1*exp(-x*x*Rg1*Rg1/3)
243        ans += B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
244        ans += G2*exp(-x*x*Rg2*Rg2/3)
245        ans += B2*(erf2^3/x)^Pow2
246        ans *= scale
247        ans += bkg
248       
249        Return(ans)
250End
251
252
253Function ThreeLevel(w,x) :FitFunc
254        Wave w
255        Variable x
256       
257        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
258        Variable G3,Rg3,B3,Pow3,erf3
259        Variable erf1,erf2,prec=1e-15,scale
260       
261        //Rsub = Rs
262        scale = w[0]
263        G1 = w[1]       //equivalent to I(0)
264        Rg1 = w[2]
265        B1 = w[3]
266        Pow1 = w[4]
267        G2 = w[5]
268        Rg2 = w[6]
269        B2 = w[7]
270        Pow2 = w[8]
271        G3 = w[9]
272        Rg3 = w[10]
273        B3 = w[11]
274        Pow3 = w[12]
275        bkg = w[13]
276       
277        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
278        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
279        erf3 = erf( (x*Rg3/sqrt(6)) ,prec)
280        //Print erf1
281       
282        ans = G1*exp(-x*x*Rg1*Rg1/3) + B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
283        ans += G2*exp(-x*x*Rg2*Rg2/3) + B2*exp(-x*x*Rg3*Rg3/3)*(erf2^3/x)^Pow2
284        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*(erf3^3/x)^Pow3
285        ans *= scale
286        ans += bkg
287       
288        Return(ans)
289End
290
291Function FourLevel(w,x) :FitFunc
292        Wave w
293        Variable x
294       
295        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
296        Variable G3,Rg3,B3,Pow3,erf3
297        Variable G4,Rg4,B4,Pow4,erf4
298        Variable erf1,erf2,prec=1e-15,scale
299       
300        //Rsub = Rs
301        scale = w[0]
302        G1 = w[1]       //equivalent to I(0)
303        Rg1 = w[2]
304        B1 = w[3]
305        Pow1 = w[4]
306        G2 = w[5]
307        Rg2 = w[6]
308        B2 = w[7]
309        Pow2 = w[8]
310        G3 = w[9]
311        Rg3 = w[10]
312        B3 = w[11]
313        Pow3 = w[12]
314        G4 = w[13]
315        Rg4 = w[14]
316        B4 = w[15]
317        Pow4 = w[16]
318        bkg = w[17]
319       
320        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
321        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
322        erf3 = erf( (x*Rg3/sqrt(6)) ,prec)
323        erf4 = erf( (x*Rg4/sqrt(6)) ,prec)
324       
325        ans = G1*exp(-x*x*Rg1*Rg1/3) + B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
326        ans += G2*exp(-x*x*Rg2*Rg2/3) + B2*exp(-x*x*Rg3*Rg3/3)*(erf2^3/x)^Pow2
327        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*exp(-x*x*Rg4*Rg4/3)*(erf3^3/x)^Pow3
328        ans += G4*exp(-x*x*Rg4*Rg4/3) + B4*(erf4^3/x)^Pow4
329        ans *= scale
330        ans += bkg
331       
332        Return(ans)
333End
334
335Function SmearedOneLevel(w,x) :FitFunc
336        Wave w
337        Variable x
338       
339//      Variable timer=StartMSTimer
340        Variable ans
341        SVAR sq = gSig_Q
342        SVAR qb = gQ_bar
343        SVAR sh = gShadow
344        SVAR gQ = gQVals
345       
346        ans = Smear_Model_20(OneLevel,$sq,$qb,$sh,$gQ,w,x)     
347
348//      Print "HS elapsed time(s) = ",StopMSTimer(timer)*1e-6
349        return(ans)
350End
351
352Function SmearedTwoLevel(w,x) :FitFunc
353        Wave w
354        Variable x
355       
356        Variable ans
357        SVAR sq = gSig_Q
358        SVAR qb = gQ_bar
359        SVAR sh = gShadow
360        SVAR gQ = gQVals
361       
362        ans = Smear_Model_20(TwoLevel,$sq,$qb,$sh,$gQ,w,x)     
363
364        return(ans)
365End
366
367Function SmearedThreeLevel(w,x) :FitFunc
368        Wave w
369        Variable x
370       
371        Variable ans
372        SVAR sq = gSig_Q
373        SVAR qb = gQ_bar
374        SVAR sh = gShadow
375        SVAR gQ = gQVals
376       
377        ans = Smear_Model_20(ThreeLevel,$sq,$qb,$sh,$gQ,w,x)   
378
379        return(ans)
380End
381
382Function SmearedFourLevel(w,x) :FitFunc
383        Wave w
384        Variable x
385       
386        Variable ans
387        SVAR sq = gSig_Q
388        SVAR qb = gQ_bar
389        SVAR sh = gShadow
390        SVAR gQ = gQVals
391       
392        ans = Smear_Model_20(FourLevel,$sq,$qb,$sh,$gQ,w,x)     
393
394        return(ans)
395End
Note: See TracBrowser for help on using the repository browser.