source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/CSParallelepiped_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: 9.2 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       
53        AddModelToStrings("CSParallpiped","coef_CSParallpiped","parameters_CSParallpiped","CSParallpiped")
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:
95        AddModelToStrings("SmearedCSParallpiped","smear_coef_CSParallpiped","smear_parameters_CSParallpiped","CSParallpiped")
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.