source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Parallelepiped_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: 6.7 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 (A^1) for model: "
33        Prompt qmax "Enter maximum q-value (A^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 (A)","B (A)","Longest Edge C (A)","SLD particle (A^-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 (A\\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","parameters_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 (A)","B (A)","Longest Edge C (A)","SLD particle (A^-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 (A\\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","smear_parameters_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.