|
tp@0
|
1 function Q = EDB1quadstep(a,b,fa,fc,fb,tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec)
|
|
tp@0
|
2 % EDB1quadstep - Recursive numerical integration routine for EDB1wedge1st_int.
|
|
tp@0
|
3 % Modified from Matlabs QUADSTEP function. The EDB1quadstep version
|
|
tp@0
|
4 % solves all time steps at the same time.
|
|
tp@0
|
5 %
|
|
tp@0
|
6 % Uses no special function
|
|
tp@0
|
7 %
|
|
tp@0
|
8 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
9 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
10 %
|
|
tp@0
|
11 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
12 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
13 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
14 %
|
|
tp@0
|
15 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
16 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
17 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
18 %
|
|
tp@0
|
19 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
20 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
21 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
22 % Peter Svensson (svensson@iet.ntnu.no) 20040321
|
|
tp@0
|
23 %
|
|
tp@0
|
24 % EDB1quadstep(a,b,fa,fc,fb,tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec)
|
|
tp@0
|
25
|
|
tp@0
|
26 % Evaluate integrand twice in interior of subinterval [a,b].
|
|
tp@0
|
27 c = (a + b)/2;
|
|
tp@0
|
28
|
|
tp@0
|
29 nslots = length(a);
|
|
tp@0
|
30
|
|
tp@0
|
31 % We use the matrix x for both the x-interval values and the corresponding
|
|
tp@0
|
32 % y-values
|
|
tp@0
|
33
|
|
tp@0
|
34 x = [(a + c)/2, (c + b)/2];
|
|
tp@0
|
35 x = reshape(x,2*nslots,1);
|
|
tp@0
|
36
|
|
tp@0
|
37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
tp@0
|
38 %% Instead of using a function call, we plug in the entire function here!
|
|
tp@0
|
39 %% x = betaoverml_all(x,rs,rr,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
40
|
|
tp@0
|
41 ml = sqrt( (x-zs).^2 + rs^2 ).*sqrt( (x-zr).^2 + rr^2 );
|
|
tp@0
|
42 %------------------------------------------------
|
|
tp@0
|
43 % We would have liked to use the following code:
|
|
tp@0
|
44 % y = (ml + (zvec-zs).*(zvec-zr))/(rs*rr);
|
|
tp@0
|
45 % expnyeta = (sqrt(y.^2-1) + y).^ny;
|
|
tp@0
|
46 % coshnyeta = 0.5*( expnyeta + 1./expnyeta);
|
|
tp@0
|
47 % but to save three vectors, zvec will be re-used for
|
|
tp@0
|
48 % y and expnyeta and coshnyeta
|
|
tp@0
|
49
|
|
tp@0
|
50 zvec = (ml + (x-zs).*(x-zr))/(rs*rr);
|
|
tp@0
|
51 zvec = (real(sqrt(zvec.^2-1)) + zvec).^ny;
|
|
tp@0
|
52 zvec = 0.5*( zvec + 1./zvec);
|
|
tp@0
|
53 x = sinnyfivec(1)./(zvec - cosnyfivec(1));
|
|
tp@0
|
54 x = x + sinnyfivec(2)./(zvec - cosnyfivec(2));
|
|
tp@0
|
55 x = x + sinnyfivec(3)./(zvec - cosnyfivec(3));
|
|
tp@0
|
56 x = x + sinnyfivec(4)./(zvec - cosnyfivec(4));
|
|
tp@0
|
57
|
|
tp@0
|
58 x = x./ml;
|
|
tp@0
|
59
|
|
tp@0
|
60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
tp@0
|
61
|
|
tp@0
|
62 x = reshape(x,nslots,2);
|
|
tp@0
|
63
|
|
tp@0
|
64 fd = x(:,1);
|
|
tp@0
|
65 fe = x(:,2);
|
|
tp@0
|
66 ftemp = fa + fb + 2*fc;
|
|
tp@0
|
67
|
|
tp@0
|
68 %-----------------------------------------
|
|
tp@0
|
69 % We use the output vector Q for another purpose
|
|
tp@0
|
70
|
|
tp@0
|
71 tempvec = (b-a)/6;
|
|
tp@0
|
72
|
|
tp@0
|
73 % Three point Simpson's rule.
|
|
tp@0
|
74 Q1 = tempvec.*(ftemp + 2*fc);
|
|
tp@0
|
75
|
|
tp@0
|
76 % Five point double Simpson's rule.
|
|
tp@0
|
77 %Q2 = (h/12).*(fa + 4*fd + 2*fc + 4*fe + fb);
|
|
tp@0
|
78
|
|
tp@0
|
79 % We would have liked to use the following code:
|
|
tp@0
|
80 % % % Q2 = tempvec/2.*(ftemp + 4*fd + 4*fe);
|
|
tp@0
|
81 % % % % One step of Romberg extrapolation.
|
|
tp@0
|
82 % % % Q = Q2 + (Q2 - Q1)/15;
|
|
tp@0
|
83 % ...... but we save one vector by re-using ftemp
|
|
tp@0
|
84 % for Q2
|
|
tp@0
|
85
|
|
tp@0
|
86 ftemp = tempvec/2.*(ftemp + 4*fd + 4*fe);
|
|
tp@0
|
87 % ftemp = tempvec/2.*(ftemp + 4*x(:,1) + 4*x(:,2));
|
|
tp@0
|
88 % One step of Romberg extrapolation.
|
|
tp@0
|
89 Q = ftemp + (ftemp - Q1)/15;
|
|
tp@0
|
90
|
|
tp@0
|
91 %-----------------------------------------
|
|
tp@0
|
92 % Check if all or some of the sample intervals
|
|
tp@0
|
93 % have reached the tolerance.
|
|
tp@0
|
94
|
|
tp@0
|
95 tempvec = find(abs(ftemp-Q)>tol);
|
|
tp@0
|
96
|
|
tp@0
|
97 % Check accuracy of integral over this subinterval.
|
|
tp@0
|
98
|
|
tp@0
|
99 if isempty(tempvec)
|
|
tp@0
|
100 return
|
|
tp@0
|
101
|
|
tp@0
|
102 % Subdivide into two subintervals.
|
|
tp@0
|
103 else
|
|
tp@0
|
104 Qac = EDB1quadstep(a(tempvec),c(tempvec),fa(tempvec),fd(tempvec),fc(tempvec),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
105 Qcb = EDB1quadstep(c(tempvec),b(tempvec),fc(tempvec),fe(tempvec),fb(tempvec),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
106 Q(tempvec) = Qac + Qcb;
|
|
tp@0
|
107 end
|
