# source:sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/CSParallelepiped_v40.ipf@504

Last change on this file since 504 was 504, checked in by ajj, 14 years ago

File size: 9.1 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// Core Shell version DS 1stJan/7thJan 2008: accounts for contribution from rims on sides A and B of diff SLDs
24//Shell on longest edge C is not included in this version
25//
26//Modified for IgorVersion 6.0 - ACH 1/7/09
27////////////////////////////////////////////////////
28
29//this macro sets up all the necessary parameters and waves that are
30//needed to calculate the model function.
31//
32Proc PlotCSParallpiped(num,qmin,qmax)
33        Variable num=100, qmin=.001, qmax=.7
34        Prompt num "Enter number of data points for model: "
35        Prompt qmin "Enter minimum q-value (A^1) for model: "
36        Prompt qmax "Enter maximum q-value (A^1) for model: "
37        //
38        Make/O/D/n=(num) xwave_CSParallpiped, ywave_CSParallpiped
39        xwave_CSParallpiped =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
40        Make/O/D coef_CSParallpiped = {1,35,75,400,10,10,10,2e-6,4e-6,2e-6,-1e-6,6e-6,0.06}                     //CH#2
41        make/o/t parameters_CSParallpiped = {"Scale Factor","Shortest Edge A (A)","B (A)","Longest Edge C (A)","Rim A ()", "Rim B ()","Rim C ()", "SLD A(A^-2)", "SLD B(A^-2)", "SLD C(A^-2)", "SLD P(A^-2)", "SLD Solv(A^-2)", "Incoherent Bgd (cm-1)"}        //CH#3
42        Edit parameters_CSParallpiped, coef_CSParallpiped
43
44        Variable/G root:g_CSParallpiped
45        g_CSParallpiped := CSParallpiped(coef_CSParallpiped,ywave_CSParallpiped, xwave_CSParallpiped)
46        Display ywave_CSParallpiped vs xwave_CSParallpiped
47        ModifyGraph marker=29, msize=2, mode=4
48        ModifyGraph log=1
49        Label bottom "q (A\\S-1\\M) "
50        Label left "I(q) (cm\\S-1\\M)"
51        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
52
54//
55End
56
57//
58//this macro sets up all the necessary parameters and waves that are
59//needed to calculate the  smeared model function.
60//
61//no input parameters are necessary, it MUST use the experimental q-values
62// from the experimental data read in from an AVE/QSIG data file
63////////////////////////////////////////////////////
64Proc PlotSmearedCSParallpiped(str)
65        String str
66        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
67
68        // if no gQvals wave, data must not have been loaded => abort
69        if(ResolutionWavesMissingDF(str))
70                Abort
71        endif
72
73        SetDataFolder \$("root:"+str)
74
75        // Setup parameter table for model function
76        Make/O/D smear_coef_CSParallpiped = {1,35,75,400,10,10,10,2e-6,4e-6,2e-6,-1e-6,6e-6,0.06}               //CH#4
77        make/o/t smear_parameters_CSParallpiped = {"Scale Factor","Shortest Edge A (A)","B (A)","Longest Edge C (A)","Rim A ()", "Rim B ()","Rim C ()", "SLD A(A^-2)", "SLD B(A^-2)", "SLD C(A^-2)", "SLD P(A^-2)", "SLD Solv(A^-2)", "Incoherent Bgd (cm-1)"}
78        Edit smear_parameters_CSParallpiped,smear_coef_CSParallpiped                                    //display parameters in a table
79
80        // output smeared intensity wave, dimensions are identical to experimental QSIG values
81        // make extra copy of experimental q-values for easy plotting
82        Duplicate/O \$(str+"_q") smeared_CSParallpiped,smeared_qvals                             //
83        SetScale d,0,0,"1/cm",smeared_CSParallpiped                                                     //
84
85        Variable/G gs_CSParallpiped=0
86        gs_CSParallpiped := fSmearedCSParallpiped(smear_coef_CSParallpiped,smeared_CSParallpiped,smeared_qvals) //this wrapper fills the STRUCT
87
88        Display smeared_CSParallpiped vs smeared_qvals                                                                  //
89        ModifyGraph log=1,marker=29,msize=2,mode=4
90        Label bottom "q (A\\S-1\\M)"
91        Label left "I(q) (cm\\S-1\\M)"
92        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
93
94        SetDataFolder root:
96End
97
98Function CSParallpiped(cw,yw,xw) : FitFunc
99        Wave cw,yw,xw
100
101#if exists("CSParallpipedX")
102        yw = CSParallpipedX(cw,xw)
103#else
104        yw = fCSParallpiped(cw,xw)
105#endif
106        return(0)
107End
108
109
110// calculates the form factor of a rectangular solid
111// - a double integral - choose points wisely
112//
113Function fCSParallpiped(w,x) : FitFunc
114        Wave w
115        Variable x
116//       Input (fitting) variables are:
117        //[0] scale factor
118        //[1] Edge A (A)
119        //[2] Edge B (A)
120        //[3] Edge C (A)
121        //[4] Rim A (A)
122        //[5] Rim B (A)
123        //[6] Rim C(A)
124        //[7] Rim A SLD(A^-2)
125        //[8] Rim B SLD (A^-2)
126        //[9] Rim C SLD (A^-2)  not included at the moment
127        //[10] PPcore SLD (A^-2)
128        //[11] Solvent SLD (A^-2)
129        //[12] incoherent background (cm^-1)
130//      give them nice names
131        Variable scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg,inten,qq,ii,arg,mu
132        scale = w[0]
133        aa = w[1]
134        bb = w[2]
135        cc = w[3]
136        ta  = w[4]
137        tb  = w[5]
138        tc  = w[6]   // is 0 at the moment
139        rhoA=w[7]   //rim A SLD
140        rhoB=w[8]   //rim B SLD
141        rhoC=w[9]    //rim C SLD
142        rhoP = w[10]   //Parallelpiped core SLD
143       rhosolv=w[11]  // Solvent SLD
144         bkg = w[12]
145
146//      mu = bb*x               //scale in terms of B
147//      aa = aa/bb
148//      cc =cc/bb
149//      ta= (aa+2*ta)/(b+2*tb)
150//      tc= (cc+2*tc)/(b+2*tb)
151
152        inten = IntegrateFn20(CSPP_Outer,0,1,w,x)
153//      inten = IntegrateFn76(PP_Outer,0,1,w,x)
154
155        inten /= (aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc)            //divide by outer volume (=Volume of core+edges)
156//      inten /= (aa*bb*cc)                                                                             // divide by volume here since FF integral has Vol^2 inside
157        inten *= 1e8            //convert to cm^-1
158//      inten *= contr*contr
159        inten *= scale
160        inten += bkg
161
162        Return (inten)
163End
164
165// outer integral
166
167// x is the q-value - remember that "mu" in the notation = B*Q
168Function CSPP_Outer(w,x,dum)
169        Wave w
170        Variable x,dum
171
172        Variable retVal,mu,mu1,aa,bb,cc,mudum, mudum1, arg
173        aa = w[1]
174        bb = w[2]
175        cc = w[3]
176//      ta  = w[4]
177//      tb  = w[5]
178//      tc  = w[6]
179
180        mu= bb*x
181//      mu1=(bb+2*tb)*x                                 //mu1 needed for including edge C contribution also: none at the moment
182        mudum = mu*sqrt(1-dum^2)
183//      mudum1 = mu1*sqrt(1-dum^2)
184        retval = IntegrateFn20(CSPP_inner,0,1,w,mudum)
185//      retval = IntegrateFn76(CSPP_inner,0,1,w,mudum)
186
187        cc = cc/bb
188        arg = mu*cc*dum/2
189        if(arg==0)
190                retval *= 1
191        else
192                retval *= (sin(arg)/arg)*(sin(arg)/arg)
193        endif
194
195        return(retVal)
196End
197
198
199//returns the integrand of the inner integral
200Function CSPP_Inner(w,mu,uu)
201        Wave w
202        Variable mu, uu    //mu1 needed for including edge C contribution also
203        Variable aa,bb,cc, ta,tb,tc, Vin,Vot,V1,V2,V3,rhoA,rhoB,rhoC, rhoP, rhosolv,dr0, drA,drB, drC,retVal,arg0,arg1,arg2,arg3,arg4,arg5,t0,t1,t2, t3, t4,t5
204                                                                                                                                                                                                        //local variables
205        //NVAR mu = root:gEvalQval              //already has been converted to S=2*pi*q
206        aa = w[1]
207        bb = w[2]
208        cc = w[3]
209         ta = w[4]
210         tb = w[5]
211         tc = w[6]
212        rhoA=w[7]
213        rhoB=w[8]
214        rhoC=w[9]
215        rhoP=w[10]
216        rhosolv=w[11]
217        dr0=rhoP-rhosolv
218        drA=rhoA-rhosolv
219        drB=rhoB-rhosolv
220        drC=rhoC-rhosolv
221         Vin=(aa*bb*cc)
222         Vot=(aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc)
223         V1=(2*ta*bb*cc)   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
224         V2=(2*aa*tb*cc)  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
225//       V3=(aa*bb*cc+2*aa*bb*tc)
226         aa = aa/bb
227//       bb = bb/bb
228         ta=(aa+2*ta)/bb
229         tb=(aa+2*tb)/bb
230
231        //Mu*(1-x^2)^(0.5)
232        //handle arg=0 separately, as sin(t)/t -> 1 as t->0
233
234//      arg0 = mu*cc*dum/2
235        arg1 = (mu*aa/2)*sin(Pi*uu/2)
236        arg2 = (mu/2)*cos(Pi*uu/2)
237        arg3=  (mu*ta/2)*sin(Pi*uu/2)
238        arg4=  (mu*tb/2)*cos(Pi*uu/2)
239
240
241//      if(arg0 ==0)
242//              t0=1
243//      else
244//              t0=sin(arg0)/arg0
245//      endif
246        if(arg1==0)
247                t1 = 1
248        else
249                t1 = (sin(arg1)/arg1)                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)
250        endif
251        if(arg2==0)
252                t2 = 1
253        else
254                t2 = (sin(arg2)/arg2)           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)
255        endif
256        if(arg3==0)
257                t3 = 1
258        else
259                t3 = sin(arg3)/arg3
260        endif
261        if(arg4==0)
262                t4 = 1
263        else
264                t4 = sin(arg4)/arg4
265        endif
266
267//      retval =( (dr0*dr0)*(t1*t1)*(t2*t2)*Vin + (drA*drA)*(t3-t1)*(t3-t1)*(t2*t2)*V1+ (drB*drB)*(t1*t1)*(t4-t2)*(t4-t2)*V2  )   // doesnot include contribution from edge C  Incorrect FF
268        retval =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )   //  correct FF : square of sum of phase factors
269//      retval =t1*t2* dr0*dr0*Vin*Vin    //*( dr0*t1*t2*Vin )   //test case of original PP with no rims
270        return(retVal);
271End
272
273
274// this is all there is to the smeared calculation!
275Function fSmearedCSParallpiped(coefW,yW,xW)
276        Wave coefW,yW,xW
277
278        String str = getWavesDataFolder(yW,0)
279        String DF="root:"+str+":"
280
281        WAVE resW = \$(DF+str+"_res")
282
283        STRUCT ResSmearAAOStruct fs
284        WAVE fs.coefW = coefW
285        WAVE fs.yW = yW
286        WAVE fs.xW = xW
287        WAVE fs.resW = resW
288
289        Variable err
290        err = SmearedCSParallpiped(fs)
291
292        return (0)
293End
294
295// this is all there is to the smeared calculation!
296Function SmearedCSParallpiped(s) :FitFunc
297        Struct ResSmearAAOStruct &s
298
299//      the name of your unsmeared model (AAO) is the first argument
300        Smear_Model_20(CSParallpiped,s.coefW,s.xW,s.yW,s.resW)
301
302        return(0)
303End
Note: See TracBrowser for help on using the repository browser.