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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

File size: 6.4 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,6e-6,0}                     //CH#2
38        make/o/t parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","Contrast (^-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//
50End
51
52// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
53Proc PlotSmearedParallelepiped(str)                                                             
54        String str
55        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
56       
57        // if any of the resolution waves are missing => abort
58        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
59                Abort
60        endif
61       
62        SetDataFolder $("root:"+str)
63       
64        // Setup parameter table for model function
65        Make/O/D smear_coef_Parallelepiped = {1,35,75,400,6e-6,0}               //CH#4
66        make/o/t smear_parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","Contrast (^-2)","Incoherent Bgd (cm-1)"}
67        Edit smear_parameters_Parallelepiped,smear_coef_Parallelepiped                                  //display parameters in a table
68       
69        // output smeared intensity wave, dimensions are identical to experimental QSIG values
70        // make extra copy of experimental q-values for easy plotting
71        Duplicate/O $(str+"_q") smeared_Parallelepiped,smeared_qvals                            //
72        SetScale d,0,0,"1/cm",smeared_Parallelepiped                                                    //
73                                       
74        Variable/G gs_Parallelepiped=0
75        gs_Parallelepiped := fSmearedParallelepiped(smear_coef_Parallelepiped,smeared_Parallelepiped,smeared_qvals)     //this wrapper fills the STRUCT
76       
77        Display smeared_Parallelepiped vs smeared_qvals                                                                 //
78        ModifyGraph log=1,marker=29,msize=2,mode=4
79        Label bottom "q (\\S-1\\M)"
80        Label left "I(q) (cm\\S-1\\M)"
81        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
82       
83        SetDataFolder root:
84End
85
86
87
88//AAO version, uses XOP if available
89// simply calls the original single point calculation with
90// a wave assignment (this will behave nicely if given point ranges)
91Function Parallelepiped(cw,yw,xw) : FitFunc
92        Wave cw,yw,xw
93       
94#if exists("ParallelepipedX")
95        yw = ParallelepipedX(cw,xw)
96#else
97        yw = fParallelepiped(cw,xw)
98#endif
99        return(0)
100End
101
102// calculates the form factor of a rectangular solid
103// - a double integral - choose points wisely
104//
105Function fParallelepiped(w,x) : FitFunc
106        Wave w
107        Variable x
108//       Input (fitting) variables are:
109        //[0] scale factor
110        //[1] Edge A (A)
111        //[2] Edge B (A)
112        //[3] Edge C (A)
113        //[4] contrast (A^-2)
114        //[5] incoherent background (cm^-1)
115//      give them nice names
116        Variable scale,aa,bb,cc,contr,bkg,inten,qq,ii,arg,mu
117        scale = w[0]
118        aa = w[1]
119        bb = w[2]
120        cc = w[3]
121        contr = w[4]
122        bkg = w[5]
123       
124//      mu = bb*x               //scale in terms of B
125//      aa = aa/bb
126//      cc = cc/bb
127       
128        inten = IntegrateFn20(PP_Outer,0,1,w,x)
129//      inten = IntegrateFn76(PP_Outer,0,1,w,x)
130       
131        inten *= aa*cc*bb               //multiply by volume
132        inten *= 1e8            //convert to cm^-1
133        inten *= contr*contr
134        inten *= scale
135        inten += bkg
136       
137        Return (inten)
138End
139
140// outer integral
141
142// x is the q-value - remember that "mu" in the notation = B*Q
143Function PP_Outer(w,x,dum)
144        Wave w
145        Variable x,dum
146       
147        Variable retVal,mu,aa,bb,cc,mudum,arg
148        aa = w[1]
149        bb = w[2]
150        cc = w[3]
151        mu = bb*x
152       
153        mudum = mu*sqrt(1-dum^2)
154        retval = IntegrateFn20(PP_inner,0,1,w,mudum)
155//      retval = IntegrateFn76(PP_inner,0,1,w,mudum)
156       
157        cc = cc/bb
158        arg = mu*cc*dum/2
159        if(arg==0)
160                retval *= 1
161        else
162                retval *= sin(arg)*sin(arg)/arg/arg
163        endif
164       
165        return(retVal)
166End
167
168//returns the integrand of the inner integral
169Function PP_Inner(w,mu,uu)
170        Wave w
171        Variable mu,uu
172       
173        Variable aa,bb,retVal,arg1,arg2,tmp1,tmp2
174       
175        //NVAR mu = root:gEvalQval              //already has been converted to S=2*pi*q
176        aa = w[1]
177        bb = w[2]
178        aa = aa/bb
179       
180        //Mu*(1-x^2)^(0.5)
181       
182        //handle arg=0 separately, as sin(t)/t -> 1 as t->0
183        arg1 = (mu/2)*cos(Pi*uu/2)
184        arg2 = (mu*aa/2)*sin(Pi*uu/2)
185        if(arg1==0)
186                tmp1 = 1
187        else
188                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1
189        endif
190        if(arg2==0)
191                tmp2 = 1
192        else
193                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2
194        endif
195        retval = tmp1*tmp2
196       
197        return(retVal)
198End
199
200//wrapper to calculate the smeared model as an AAO-Struct
201// fills the struct and calls the ususal function with the STRUCT parameter
202//
203// used only for the dependency, not for fitting
204//
205Function fSmearedParallelepiped(coefW,yW,xW)
206        Wave coefW,yW,xW
207       
208        String str = getWavesDataFolder(yW,0)
209        String DF="root:"+str+":"
210       
211        WAVE resW = $(DF+str+"_res")
212       
213        STRUCT ResSmearAAOStruct fs
214        WAVE fs.coefW = coefW   
215        WAVE fs.yW = yW
216        WAVE fs.xW = xW
217        WAVE fs.resW = resW
218       
219        Variable err
220        err = SmearedParallelepiped(fs)
221       
222        return (0)
223End
224
225// this is all there is to the smeared calculation!
226Function SmearedParallelepiped(s) :FitFunc
227        Struct ResSmearAAOStruct &s
228
229//      the name of your unsmeared model (AAO) is the first argument
230        Smear_Model_20(Parallelepiped,s.coefW,s.xW,s.yW,s.resW)
231
232        return(0)
233End
Note: See TracBrowser for help on using the repository browser.