source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/Parallelepiped_v40.ipf @ 345

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

Merging NewModels_2006 back into where it belongs

File size: 6.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////////
5//
6// calculates the scattering from a rectangular solid
7// i.e. a parallelepiped with sides a < b < c
8//
9// - the user must make sure that the constraints are not violated
10// otherwise the calculation will not be correct
11//
12// From: Mittelbach and Porod, Acta Phys. Austriaca 14 (1961) 185-211.
13//                              equations (1), (13), and (14) (in German!)
14//
15// note that the equations listed in Feigin and Svergun appears
16// to be wrong - they use equation (12), which does not appear to
17// be a complete orientational average (?)
18//
19// a double integral is used, both using Gaussian quadrature
20// routines that are now included with GaussUtils
21// 20-pt quadrature appears to be enough, 76 pt is available
22// by changing the function calls
23//
24////////////////////////////////////////////////////
25
26//this macro sets up all the necessary parameters and waves that are
27//needed to calculate the model function.
28//
29Proc PlotParallelepiped(num,qmin,qmax)
30        Variable num=100, qmin=.001, qmax=.7
31        Prompt num "Enter number of data points for model: "
32        Prompt qmin "Enter minimum q-value (^1) for model: "
33        Prompt qmax "Enter maximum q-value (^1) for model: "
34        //
35        Make/O/D/n=(num) xwave_Parallelepiped, ywave_Parallelepiped
36        xwave_Parallelepiped =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
37        Make/O/D coef_Parallelepiped = {1,35,75,400,1e-6,6.3e-6,0}                      //CH#2
38        make/o/t parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","SLD particle (^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}     //CH#3
39        Edit parameters_Parallelepiped, coef_Parallelepiped
40       
41        Variable/G root:g_Parallelepiped
42        g_Parallelepiped := Parallelepiped(coef_Parallelepiped,ywave_Parallelepiped, xwave_Parallelepiped)
43        Display ywave_Parallelepiped vs xwave_Parallelepiped
44        ModifyGraph marker=29, msize=2, mode=4
45        ModifyGraph log=1
46        Label bottom "q (\\S-1\\M) "
47        Label left "I(q) (cm\\S-1\\M)"
48        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
49       
50        AddModelToStrings("Parallelepiped","coef_Parallelepiped","Parallelepiped")
51//
52End
53
54// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
55Proc PlotSmearedParallelepiped(str)                                                             
56        String str
57        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
58       
59        // if any of the resolution waves are missing => abort
60        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
61                Abort
62        endif
63       
64        SetDataFolder $("root:"+str)
65       
66        // Setup parameter table for model function
67        Make/O/D smear_coef_Parallelepiped = {1,35,75,400,1e-6,6.3e-6,0}                //CH#4
68        make/o/t smear_parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","SLD particle (^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
69        Edit smear_parameters_Parallelepiped,smear_coef_Parallelepiped                                  //display parameters in a table
70       
71        // output smeared intensity wave, dimensions are identical to experimental QSIG values
72        // make extra copy of experimental q-values for easy plotting
73        Duplicate/O $(str+"_q") smeared_Parallelepiped,smeared_qvals                            //
74        SetScale d,0,0,"1/cm",smeared_Parallelepiped                                                    //
75                                       
76        Variable/G gs_Parallelepiped=0
77        gs_Parallelepiped := fSmearedParallelepiped(smear_coef_Parallelepiped,smeared_Parallelepiped,smeared_qvals)     //this wrapper fills the STRUCT
78       
79        Display smeared_Parallelepiped vs smeared_qvals                                                                 //
80        ModifyGraph log=1,marker=29,msize=2,mode=4
81        Label bottom "q (\\S-1\\M)"
82        Label left "I(q) (cm\\S-1\\M)"
83        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
84       
85        SetDataFolder root:
86        AddModelToStrings("SmearedParallelepiped","smear_coef_Parallelepiped","Parallelepiped")
87End
88
89
90
91//AAO version, uses XOP if available
92// simply calls the original single point calculation with
93// a wave assignment (this will behave nicely if given point ranges)
94Function Parallelepiped(cw,yw,xw) : FitFunc
95        Wave cw,yw,xw
96       
97#if exists("ParallelepipedX")
98        yw = ParallelepipedX(cw,xw)
99#else
100        yw = fParallelepiped(cw,xw)
101#endif
102        return(0)
103End
104
105// calculates the form factor of a rectangular solid
106// - a double integral - choose points wisely
107//
108Function fParallelepiped(w,x) : FitFunc
109        Wave w
110        Variable x
111//       Input (fitting) variables are:
112        //[0] scale factor
113        //[1] Edge A (A)
114        //[2] Edge B (A)
115        //[3] Edge C (A)
116        //[4] contrast (A^-2)
117        //[5] incoherent background (cm^-1)
118//      give them nice names
119        Variable scale,aa,bb,cc,contr,bkg,inten,qq,ii,arg,mu,sldp,slds
120        scale = w[0]
121        aa = w[1]
122        bb = w[2]
123        cc = w[3]
124        sldp = w[4]
125        slds = w[5]
126        bkg = w[6]
127       
128        contr = sldp - slds
129//      mu = bb*x               //scale in terms of B
130//      aa = aa/bb
131//      cc = cc/bb
132       
133        inten = IntegrateFn20(PP_Outer,0,1,w,x)
134//      inten = IntegrateFn76(PP_Outer,0,1,w,x)
135       
136        inten *= aa*cc*bb               //multiply by volume
137        inten *= 1e8            //convert to cm^-1
138        inten *= contr*contr
139        inten *= scale
140        inten += bkg
141       
142        Return (inten)
143End
144
145// outer integral
146
147// x is the q-value - remember that "mu" in the notation = B*Q
148Function PP_Outer(w,x,dum)
149        Wave w
150        Variable x,dum
151       
152        Variable retVal,mu,aa,bb,cc,mudum,arg
153        aa = w[1]
154        bb = w[2]
155        cc = w[3]
156        mu = bb*x
157       
158        mudum = mu*sqrt(1-dum^2)
159        retval = IntegrateFn20(PP_inner,0,1,w,mudum)
160//      retval = IntegrateFn76(PP_inner,0,1,w,mudum)
161       
162        cc = cc/bb
163        arg = mu*cc*dum/2
164        if(arg==0)
165                retval *= 1
166        else
167                retval *= sin(arg)*sin(arg)/arg/arg
168        endif
169       
170        return(retVal)
171End
172
173//returns the integrand of the inner integral
174Function PP_Inner(w,mu,uu)
175        Wave w
176        Variable mu,uu
177       
178        Variable aa,bb,retVal,arg1,arg2,tmp1,tmp2
179       
180        //NVAR mu = root:gEvalQval              //already has been converted to S=2*pi*q
181        aa = w[1]
182        bb = w[2]
183        aa = aa/bb
184       
185        //Mu*(1-x^2)^(0.5)
186       
187        //handle arg=0 separately, as sin(t)/t -> 1 as t->0
188        arg1 = (mu/2)*cos(Pi*uu/2)
189        arg2 = (mu*aa/2)*sin(Pi*uu/2)
190        if(arg1==0)
191                tmp1 = 1
192        else
193                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1
194        endif
195        if(arg2==0)
196                tmp2 = 1
197        else
198                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2
199        endif
200        retval = tmp1*tmp2
201       
202        return(retVal)
203End
204
205//wrapper to calculate the smeared model as an AAO-Struct
206// fills the struct and calls the ususal function with the STRUCT parameter
207//
208// used only for the dependency, not for fitting
209//
210Function fSmearedParallelepiped(coefW,yW,xW)
211        Wave coefW,yW,xW
212       
213        String str = getWavesDataFolder(yW,0)
214        String DF="root:"+str+":"
215       
216        WAVE resW = $(DF+str+"_res")
217       
218        STRUCT ResSmearAAOStruct fs
219        WAVE fs.coefW = coefW   
220        WAVE fs.yW = yW
221        WAVE fs.xW = xW
222        WAVE fs.resW = resW
223       
224        Variable err
225        err = SmearedParallelepiped(fs)
226       
227        return (0)
228End
229
230// this is all there is to the smeared calculation!
231Function SmearedParallelepiped(s) :FitFunc
232        Struct ResSmearAAOStruct &s
233
234//      the name of your unsmeared model (AAO) is the first argument
235        Smear_Model_20(Parallelepiped,s.coefW,s.xW,s.yW,s.resW)
236
237        return(0)
238End
Note: See TracBrowser for help on using the repository browser.