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

Last change on this file since 91 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

File size: 5.5 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////////
4//
5// calculates the scattering from a rectangular solid
6// i.e. a parallelepiped with sides a < b < c
7//
8// - the user must make sure that the constraints are not violated
9// otherwise the calculation will not be correct
10//
11// From: Mittelbach and Porod, Acta Phys. Austriaca 14 (1961) 185-211.
12//                              equations (1), (13), and (14) (in German!)
13//
14// note that the equations listed in Feigin and Svergun appears
15// to be wrong - they use equation (12), which does not appear to
16// be a complete orientational average (?)
17//
18// a double integral is used, both using Gaussian quadrature
19// routines that are now included with GaussUtils
20// 20-pt quadrature appears to be enough, 76 pt is available
21// by changing the function calls
22//
23////////////////////////////////////////////////////
24
25//this macro sets up all the necessary parameters and waves that are
26//needed to calculate the model function.
27//
28Proc Plot_Parallelepiped(num,qmin,qmax)
29        Variable num=100, qmin=.001, qmax=.7
30        Prompt num "Enter number of data points for model: "
31        Prompt qmin "Enter minimum q-value (^1) for model: "
32        Prompt qmax "Enter maximum q-value (^1) for model: "
33        //
34        Make/O/D/n=(num) xwave_Parallelepiped, ywave_Parallelepiped
35        xwave_Parallelepiped =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
36        Make/O/D coef_Parallelepiped = {1,35,75,400,6e-6,0}                     //CH#2
37        make/o/t parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","Contrast (^-2)","Incoherent Bgd (cm-1)"}      //CH#3
38        Edit parameters_Parallelepiped, coef_Parallelepiped
39        ywave_Parallelepiped  := Parallelepiped(coef_Parallelepiped, xwave_Parallelepiped)
40        Display ywave_Parallelepiped vs xwave_Parallelepiped
41        ModifyGraph marker=29, msize=2, mode=4
42        ModifyGraph log=1
43        Label bottom "q (\\S-1\\M) "
44        Label left "I(q) (cm\\S-1\\M)"
45        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
46//
47End
48
49//
50//this macro sets up all the necessary parameters and waves that are
51//needed to calculate the  smeared model function.
52//
53//no input parameters are necessary, it MUST use the experimental q-values
54// from the experimental data read in from an AVE/QSIG data file
55////////////////////////////////////////////////////
56Proc PlotSmeared_Parallelepiped()                                                               //Parallelepiped
57        // if no gQvals wave, data must not have been loaded => abort
58        if(ResolutionWavesMissing())
59                Abort
60        endif
61
62
63        // Setup parameter table for model function
64        Make/O/D smear_coef_Parallelepiped = {1,35,75,400,6e-6,0}               //CH#4
65        make/o/t smear_parameters_Parallelepiped = {"Scale Factor","Shortest Edge A ()","B ()","Longest Edge C ()","Contrast (^-2)","Incoherent Bgd (cm-1)"}
66        Edit smear_parameters_Parallelepiped,smear_coef_Parallelepiped                                  //display parameters in a table
67
68        // output smeared intensity wave, dimensions are identical to experimental QSIG values
69        // make extra copy of experimental q-values for easy plotting
70        Duplicate/O \$gQvals smeared_Parallelepiped,smeared_qvals                                //
71        SetScale d,0,0,"1/cm",smeared_Parallelepiped                                                    //
72
73        smeared_Parallelepiped := Parallelepiped_Smeared(smear_coef_Parallelepiped,\$gQvals)             // SMEARED function name
74        Display smeared_Parallelepiped vs smeared_qvals                                                                 //
75        ModifyGraph log=1,marker=29,msize=2,mode=4
76        Label bottom "q (\\S-1\\M)"
77        Label left "I(q) (cm\\S-1\\M)"
78        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
79End     // end macro
80
81// calculates the form factor of a rectangular solid
82// - a double integral - choose points wisely
83//
84Function Parallelepiped(w,x) : FitFunc
85        Wave w
86        Variable x
87//       Input (fitting) variables are:
88        //[0] scale factor
89        //[1] Edge A (A)
90        //[2] Edge B (A)
91        //[3] Edge C (A)
92        //[4] contrast (A^-2)
93        //[5] incoherent background (cm^-1)
94//      give them nice names
95        Variable scale,aa,bb,cc,contr,bkg,inten,qq,ii,arg,mu
96        scale = w[0]
97        aa = w[1]
98        bb = w[2]
99        cc = w[3]
100        contr = w[4]
101        bkg = w[5]
102
103//      mu = bb*x               //scale in terms of B
104//      aa = aa/bb
105//      cc = cc/bb
106
107        inten = IntegrateFn20(PP_Outer,0,1,w,x)
108//      inten = IntegrateFn76(PP_Outer,0,1,w,x)
109
110        inten *= aa*cc*bb               //multiply by volume
111        inten *= 1e8            //convert to cm^-1
112        inten *= contr*contr
113        inten *= scale
114        inten += bkg
115
116        Return (inten)
117End
118
119// outer integral
120
121// x is the q-value - remember that "mu" in the notation = B*Q
122Function PP_Outer(w,x,dum)
123        Wave w
124        Variable x,dum
125
126        Variable retVal,mu,aa,bb,cc,mudum,arg
127        aa = w[1]
128        bb = w[2]
129        cc = w[3]
130        mu = bb*x
131
132        mudum = mu*sqrt(1-dum^2)
133        retval = IntegrateFn20(PP_inner,0,1,w,mudum)
134//      retval = IntegrateFn76(PP_inner,0,1,w,mudum)
135
136        cc = cc/bb
137        arg = mu*cc*dum/2
138        if(arg==0)
139                retval *= 1
140        else
141                retval *= sin(arg)*sin(arg)/arg/arg
142        endif
143
144        return(retVal)
145End
146
147//returns the integrand of the inner integral
148Function PP_Inner(w,mu,uu)
149        Wave w
150        Variable mu,uu
151
152        Variable aa,bb,retVal,arg1,arg2,tmp1,tmp2
153
154        //NVAR mu = root:gEvalQval              //already has been converted to S=2*pi*q
155        aa = w[1]
156        bb = w[2]
157        aa = aa/bb
158
159        //Mu*(1-x^2)^(0.5)
160
161        //handle arg=0 separately, as sin(t)/t -> 1 as t->0
162        arg1 = (mu/2)*cos(Pi*uu/2)
163        arg2 = (mu*aa/2)*sin(Pi*uu/2)
164        if(arg1==0)
165                tmp1 = 1
166        else
167                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1
168        endif
169        if(arg2==0)
170                tmp2 = 1
171        else
172                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2
173        endif
174        retval = tmp1*tmp2
175
176        return(retVal)
177End
178
179
180// this is all there is to the smeared calculation!
181Function Parallelepiped_Smeared(w,x) :FitFunc
182        Wave w
183        Variable x
184
185        Variable ans
186        SVAR sq = gSig_Q
187        SVAR qb = gQ_bar