|
tp@0
|
1 function [ir,ninit] = EDB1wedge2nd(cylS,cylR,cylE2_r1,cylE1_r2,...
|
|
tp@0
|
2 nyveclist,edgelengthlist,dzvec,method,pathalongplane,R_irstart,bc,cair,fs)
|
|
tp@0
|
3 % EDB1wedge2nd - Gives the 2nd-order diffraction IR.
|
|
tp@0
|
4 %
|
|
tp@0
|
5 % NB!!! Only edges that are in-plane with each other can be handled with
|
|
tp@0
|
6 % this version. A small extension is needed to handle non-in-plane edge
|
|
tp@0
|
7 % pairs.
|
|
tp@0
|
8 %
|
|
tp@0
|
9 % NB!!! Known limitation/error: geometries like reflector clouds with
|
|
tp@0
|
10 % several surfaces in the same plane will give some errors: obstruction
|
|
tp@0
|
11 % checks are not done properly for edge-to-edge pairs if they are in the
|
|
tp@0
|
12 % same plane. Some removal of non-allowed edge-to-edge pairs is done: if a
|
|
tp@0
|
13 % path from one edge to another edges starts by crossing the first edge's
|
|
tp@0
|
14 % own plane, then that edge-to-edge combination is shut off. This could be
|
|
tp@0
|
15 % expressed as if a second edge is "behind" the first edge. It corresponds
|
|
tp@0
|
16 % somewhat to the check for the image source method if an image source is
|
|
tp@0
|
17 % behind or in front of its last potential reflection surface - a
|
|
tp@0
|
18 % "validity" check. But, the subsequent "obstruction" check is not done
|
|
tp@0
|
19 % yet.
|
|
tp@0
|
20 %
|
|
tp@0
|
21 % ir The 2nd order diffraction impulse response for two
|
|
tp@0
|
22 % edges that share one plane.
|
|
tp@0
|
23 % ninit The initial delay that has been cut away before the
|
|
tp@0
|
24 % beginning of ir. Note that the time zero of the ir (including ninit)
|
|
tp@0
|
25 % corresponds to R_irstart.
|
|
tp@0
|
26 %
|
|
tp@0
|
27 % with input parameters
|
|
tp@0
|
28 % cylS vector containing [rS,thetaS,zS] = cyl. coordinates of the source rel. to the first edge
|
|
tp@0
|
29 % cylR vector containing [rR,thetaR,zR] = cyl. coordinates of the receiver rel. to the second edge
|
|
tp@0
|
30 % cylE2_r1 matrix containing [rstart_r1,thetastart_r1,zstart_r1;rend_r1,thetaend_r1,zend_r1]
|
|
tp@0
|
31 % = cyl. coordinates of the start and end points of the second edge rel. to the first edge
|
|
tp@0
|
32 % cylE1_r2 matrix containing [rstart_r2,thetastart_r2,zstart_r2;rend_r2,thetaend_r2,zend_r2]
|
|
tp@0
|
33 % = cyl. coordinates of the start and end points of the first edge rel. to the second edge
|
|
tp@0
|
34 % nyveclist
|
|
tp@0
|
35 % vector with the two wedge indices of the two edges
|
|
tp@0
|
36 % edgelengthlist
|
|
tp@0
|
37 % Matrix, [2 2], with the start and end values of the two edges:
|
|
tp@0
|
38 % Row 1 has the start and end values of edge 1
|
|
tp@0
|
39 % Row 2 has the start and end values of edge 2
|
|
tp@0
|
40 % For a fully visible/active edge, the start value = 0 and the
|
|
tp@0
|
41 % end value equals the length (in meters) of the edge.
|
|
tp@0
|
42 % elemsize_2nd
|
|
tp@0
|
43 % the number of elements per wavelength at the sampling frequency.
|
|
tp@0
|
44 % Recommended value 0.5, and this value is chosen if no other value
|
|
tp@0
|
45 % has been set. One can experiment with lower values (0.25, 0.1, ...)
|
|
tp@0
|
46 % for comparisons - the error increases with frequency but if first
|
|
tp@0
|
47 % order diffraction and specular response is masking the second order
|
|
tp@0
|
48 % diffraction, low values of elemsize_2nd might be perfectly acceptable.
|
|
tp@0
|
49 % method 'n' for the New method or 'v' for Vanderkooys method
|
|
tp@0
|
50 % pathalongplane 1 or 0 depending on if the path runs along a plane
|
|
tp@0
|
51 % or not
|
|
tp@0
|
52 % R_irstart (optional)
|
|
tp@0
|
53 % a distance to which the ir time zero will correspond
|
|
tp@0
|
54 % bc (optional)
|
|
tp@0
|
55 % a matrix containing [refl11 refl12;refl21 refl22] where the values are +1 or -1
|
|
tp@0
|
56 % indicating the reflection coefficients of the four wedge planes involved.
|
|
tp@0
|
57 % +1 means rigid and -1 means soft (pressure-release)
|
|
tp@0
|
58 % refl11 is edge 1, plane 1 (refplane)
|
|
tp@0
|
59 % refl12 is edge 1, plane 2
|
|
tp@0
|
60 % refl21 is edge 2, plane 1 (refplane)
|
|
tp@0
|
61 % refl22 is edge 2, plane 2
|
|
tp@0
|
62 % default is all-rigid
|
|
tp@0
|
63 % CAIR (global)
|
|
tp@0
|
64 % the speed of sound
|
|
tp@0
|
65 % FSAMP (global)
|
|
tp@0
|
66 % the sampling frequency
|
|
tp@0
|
67 % BIGEDGESTEPMATRIX (global)
|
|
tp@0
|
68 % a matrix, size [N,2], where N = nedge1*nedge2. nedge1 is the number of edge
|
|
tp@0
|
69 % elements that edge1 is sub-divided into and nedge2 is the
|
|
tp@0
|
70 % number of edge elements that edge2 is sub-divided into. The
|
|
tp@0
|
71 % matrix contains values between 0 and 1, representing all
|
|
tp@0
|
72 % combinations of relative positions along the two edges for the
|
|
tp@0
|
73 % mid-points of these edge elements. They are organized like
|
|
tp@0
|
74 % below, if nedge1 = 10 and nedge2 = 25:
|
|
tp@0
|
75 % BIGEDGESTEPMATRIX =
|
|
tp@0
|
76 % [ 0.05 0.02;
|
|
tp@0
|
77 % 0.05 0.06;
|
|
tp@0
|
78 % 0.05 0.10;
|
|
tp@0
|
79 % ...
|
|
tp@0
|
80 % 0.05 0.98;
|
|
tp@0
|
81 % 0.15 0.02;
|
|
tp@0
|
82 % ...
|
|
tp@0
|
83 % 0.95 0.98];
|
|
tp@0
|
84 % These values are multiplied with the lengths of the respective
|
|
tp@0
|
85 % edges to give the actual position, in meters along each edge.
|
|
tp@0
|
86 %
|
|
tp@0
|
87 % Uses no special functions
|
|
tp@0
|
88 %
|
|
tp@0
|
89 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
90 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
91 %
|
|
tp@0
|
92 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
93 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
94 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
95 %
|
|
tp@0
|
96 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
97 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
98 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
99 %
|
|
tp@0
|
100 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
101 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
102 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
103 % Peter Svensson (svensson@iet.ntnu.no) 20080508
|
|
tp@0
|
104 %
|
|
tp@0
|
105 % [ir,ninit] = EDB1wedge2ndnew(cylS,cylR,cylE2_r1,cylE1_r2,...
|
|
tp@0
|
106 % nyveclist,edgelengthlist,dzvec,method,pathalongplane,BigB,R_irstart,bc,cair,fs);
|
|
tp@0
|
107
|
|
tp@0
|
108 global BIGEDGESTEPMATRIX
|
|
tp@0
|
109
|
|
tp@0
|
110 multfac = 1/(pathalongplane + 1);
|
|
tp@0
|
111
|
|
tp@0
|
112 %-----------------------------------------------------------------------
|
|
tp@0
|
113 % Extract the individual parameters
|
|
tp@0
|
114
|
|
tp@0
|
115 rS = cylS(1);
|
|
tp@0
|
116 thetaS = cylS(2);
|
|
tp@0
|
117 zS = cylS(3);
|
|
tp@0
|
118
|
|
tp@0
|
119 rR = cylR(1);
|
|
tp@0
|
120 thetaR = cylR(2);
|
|
tp@0
|
121 zR = cylR(3);
|
|
tp@0
|
122
|
|
tp@0
|
123 rE1_r1 = cylE2_r1(1,1);
|
|
tp@0
|
124 thetaE1_r1 = cylE2_r1(1,2);
|
|
tp@0
|
125 zE1_r1 = cylE2_r1(1,3);
|
|
tp@0
|
126 rE2_r1 = cylE2_r1(2,1);
|
|
tp@0
|
127 thetaE2_r1 = cylE2_r1(2,2);
|
|
tp@0
|
128 zE2_r1 = cylE2_r1(2,3);
|
|
tp@0
|
129
|
|
tp@0
|
130 rE1_r2 = cylE1_r2(1,1);
|
|
tp@0
|
131 thetaE1_r2 = cylE1_r2(1,2);
|
|
tp@0
|
132 zE1_r2 = cylE1_r2(1,3);
|
|
tp@0
|
133 rE2_r2 = cylE1_r2(2,1);
|
|
tp@0
|
134 thetaE2_r2 = cylE1_r2(2,2);
|
|
tp@0
|
135 zE2_r2 = cylE1_r2(2,3);
|
|
tp@0
|
136
|
|
tp@0
|
137 edgesareinplane = ( abs(thetaE1_r1 - thetaE2_r1) < 1e-6 );
|
|
tp@0
|
138
|
|
tp@0
|
139 swapbigmatrix = 0;
|
|
tp@0
|
140 zmid1 = (edgelengthlist(1,1)+edgelengthlist(1,2))/2;
|
|
tp@0
|
141 zmid2_re1 = ( zE1_r1 + zE2_r1 )/2;
|
|
tp@0
|
142 midpointdist = sqrt( rE1_r1.^2 + (zmid2_re1-zmid1).^2);
|
|
tp@0
|
143 reledgetoedgedist = midpointdist/max(dzvec);
|
|
tp@0
|
144
|
|
tp@0
|
145 % % % % % % % % % % % if reledgetoedgedist < 12
|
|
tp@0
|
146 % % % % % % % % % % % refinefactor = ceil(12/reledgetoedgedist);
|
|
tp@0
|
147 % % % % % % % % % % % % disp([' Refining by factor ',int2str(refinefactor)])
|
|
tp@0
|
148 % % % % % % % % % % % newndivvec = round(edgelengthlist(:,2)./dzvec.').'*refinefactor;
|
|
tp@0
|
149 % % % % % % % % % % % ivmatrix = EDB1creindexmatrix(newndivvec);
|
|
tp@0
|
150 % % % % % % % % % % % [nedgeelcombs,slask] = size(ivmatrix);
|
|
tp@0
|
151 % % % % % % % % % % % BIGE = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),:);
|
|
tp@0
|
152 % % % % % % % % % % % swapbigmatrix = 1;
|
|
tp@0
|
153 % % % % % % % % % % % dzvec = dzvec/refinefactor;
|
|
tp@0
|
154 % % % % % % % % % % % end
|
|
tp@0
|
155
|
|
tp@0
|
156 % E2E, first part
|
|
tp@0
|
157
|
|
tp@0
|
158 if swapbigmatrix == 0
|
|
tp@0
|
159 B4 = edgelengthlist(1,1) + (edgelengthlist(1,2)-edgelengthlist(1,1))*BIGEDGESTEPMATRIX(:,1); % B4 is zedge1_re1 which goes, per def., from 0 to edgelength1
|
|
tp@0
|
160 zedge2_re1 = zE1_r1 + BIGEDGESTEPMATRIX(:,2)*( zE2_r1 - zE1_r1 ); % zedge2_re1 are the z-values of edge2, expressed in edge 1
|
|
tp@0
|
161 B5 = edgelengthlist(2,1) + (edgelengthlist(2,2)-edgelengthlist(2,1))*BIGEDGESTEPMATRIX(:,2); % B5 is zedge2_re2 which goes, per def., from 0 to edgelength2
|
|
tp@0
|
162 zedge1_re2 = zE1_r2 + BIGEDGESTEPMATRIX(:,1)*( zE2_r2 - zE1_r2 );
|
|
tp@0
|
163 else
|
|
tp@0
|
164 B4 = edgelengthlist(1,1) + (edgelengthlist(1,2)-edgelengthlist(1,1))*BIGE(:,1); % B4 is zedge1_re1 which goes, per def., from 0 to edgelength1
|
|
tp@0
|
165 zedge2_re1 = zE1_r1 + BIGE(:,2)*( zE2_r1 - zE1_r1 ); % zedge2_re1 are the z-values of edge2, expressed in edge 1
|
|
tp@0
|
166 B5 = edgelengthlist(2,1) + (edgelengthlist(2,2)-edgelengthlist(2,1))*BIGE(:,2); % B5 is zedge2_re2 which goes, per def., from 0 to edgelength2
|
|
tp@0
|
167 zedge1_re2 = zE1_r2 + BIGE(:,1)*( zE2_r2 - zE1_r2 );
|
|
tp@0
|
168 end
|
|
tp@0
|
169
|
|
tp@0
|
170 % S2E
|
|
tp@0
|
171 S2Edist = sqrt( rS^2 + ( B4 - zS ).^2);
|
|
tp@0
|
172
|
|
tp@0
|
173 % E2R
|
|
tp@0
|
174 E2Rdist = sqrt( rR^2 + ( B5 - zR ).^2);
|
|
tp@0
|
175
|
|
tp@0
|
176 % E2E, second part
|
|
tp@0
|
177 if swapbigmatrix == 0
|
|
tp@0
|
178 if edgesareinplane == 0,
|
|
tp@0
|
179
|
|
tp@0
|
180 % error('ERROR: r and theta calculation for non-in-plane edges not implemented yet!')
|
|
tp@0
|
181
|
|
tp@0
|
182 % First we need to reconvert the cylindrical coordinates of
|
|
tp@0
|
183 % edge2-re1 into cartesian! Then we can use the EDB1coordtrans.
|
|
tp@0
|
184 % We define our own cartesian coord syst such that the reference
|
|
tp@0
|
185 % edge has its starting point in [0 0 0], and the x-axis is along
|
|
tp@0
|
186 % the reference plane.
|
|
tp@0
|
187
|
|
tp@0
|
188 xE1_re1 = rE1_r1*cos(thetaE1_r1);
|
|
tp@0
|
189 yE1_re1 = rE1_r1*sin(thetaE1_r1);
|
|
tp@0
|
190 xE2_re1 = rE2_r1*cos(thetaE2_r1);
|
|
tp@0
|
191 yE2_re1 = rE2_r1*sin(thetaE2_r1);
|
|
tp@0
|
192
|
|
tp@0
|
193 xvec_re1 = xE1_re1 + (xE2_re1 - xE1_re1)*BIGEDGESTEPMATRIX(:,2);
|
|
tp@0
|
194 yvec_re1 = yE1_re1 + (yE2_re1 - yE1_re1)*BIGEDGESTEPMATRIX(:,2);
|
|
tp@0
|
195 zvec_re1 = zE1_r1 + (zE2_r1 - zE1_r1)*BIGEDGESTEPMATRIX(:,2);
|
|
tp@0
|
196
|
|
tp@0
|
197 % NB!!! Below we need to check if really the full lengths of both
|
|
tp@0
|
198 % edges should be used. If S or R see only part of the respective
|
|
tp@0
|
199 % edge then the edge-to-edge contribution should be "windowed" too.
|
|
tp@0
|
200 [redge2_re1,thetaedge2_re1,znewedge2_re1] = EDB1coordtrans1([xvec_re1 yvec_re1 zvec_re1],[0 0 0;0 0 edgelengthlist(1,2)],[0 1 0]);
|
|
tp@0
|
201
|
|
tp@0
|
202 xE1_re2 = rE1_r2*cos(thetaE1_r2);
|
|
tp@0
|
203 yE1_re2 = rE1_r2*sin(thetaE1_r2);
|
|
tp@0
|
204 xE2_re2 = rE2_r2*cos(thetaE2_r2);
|
|
tp@0
|
205 yE2_re2 = rE2_r2*sin(thetaE2_r2);
|
|
tp@0
|
206
|
|
tp@0
|
207 xvec_re2 = xE1_re2 + (xE2_re2 - xE1_re2)*BIGEDGESTEPMATRIX(:,1);
|
|
tp@0
|
208 yvec_re2 = yE1_re2 + (yE2_re2 - yE1_re2)*BIGEDGESTEPMATRIX(:,1);
|
|
tp@0
|
209 zvec_re2 = zE1_r2 + (zE2_r2 - zE1_r2)*BIGEDGESTEPMATRIX(:,1);
|
|
tp@0
|
210
|
|
tp@0
|
211 % NB!!! Below we need to check if really the full lengths of both
|
|
tp@0
|
212 % edges should be used. If S or R see only part of the respective
|
|
tp@0
|
213 % edge then the edge-to-edge contribution should be "windowed" too.
|
|
tp@0
|
214 [redge1_re2,thetaedge1_re2,znewedge1_re2] = EDB1coordtrans1([xvec_re2 yvec_re2 zvec_re2],[0 0 0;0 0 edgelengthlist(2,2)],[0 1 0]);
|
|
tp@0
|
215
|
|
tp@0
|
216 else
|
|
tp@0
|
217 redge2_re1 = rE1_r1 + BIGEDGESTEPMATRIX(:,2)*( rE2_r1 - rE1_r1 );
|
|
tp@0
|
218 redge1_re2 = rE1_r2 + BIGEDGESTEPMATRIX(:,1)*( rE2_r2 - rE1_r2 );
|
|
tp@0
|
219 end
|
|
tp@0
|
220 else
|
|
tp@0
|
221
|
|
tp@0
|
222 if edgesareinplane == 0
|
|
tp@0
|
223 error('ERROR: r and theta calculation for non-in-plane edges not implemented yet, for swapbigmatrix=1!')
|
|
tp@0
|
224
|
|
tp@0
|
225 else
|
|
tp@0
|
226 redge2_re1 = rE1_r1 + BIGE(:,2)*( rE2_r1 - rE1_r1 );
|
|
tp@0
|
227 redge1_re2 = rE1_r2 + BIGE(:,1)*( rE2_r2 - rE1_r2 );
|
|
tp@0
|
228 end
|
|
tp@0
|
229 end
|
|
tp@0
|
230 if edgesareinplane == 0
|
|
tp@0
|
231 % PS 20120415: CHECK - is it correct or not for non-in-plane edges?
|
|
tp@0
|
232 % % % error('ERROR: r and theta calculation for non-in-plane edges not implemented yet!')
|
|
tp@0
|
233 % % % % We need to implement something like:
|
|
tp@0
|
234 % % % % theta = acos(r./E2Edist);
|
|
tp@0
|
235 % % % thetaedge1_re2 = thetaE1_r2 + BIGEDGESTEPMATRIX(:,1)*( thetaE2_r2 - thetaE1_r2 );
|
|
tp@0
|
236 % % % thetaedge2_re1 = thetaE1_r1 + BIGEDGESTEPMATRIX(:,2)*( thetaE2_r1 - thetaE1_r1 );
|
|
tp@0
|
237 else
|
|
tp@0
|
238 thetaedge1_re2 = thetaE1_r2;
|
|
tp@0
|
239 thetaedge2_re1 = thetaE1_r1;
|
|
tp@0
|
240 end
|
|
tp@0
|
241
|
|
tp@0
|
242 if method == 'n'
|
|
tp@0
|
243
|
|
tp@0
|
244 %-----------------------------------------------------------------------
|
|
tp@0
|
245 % In directivity function 2, DF2, we need the quantity
|
|
tp@0
|
246 %
|
|
tp@0
|
247 % ( ( ze2_re2 - ze1_re2 ) * ( ze2_re2 - zr_re2 ) + n*l )/re1_re2/rr
|
|
tp@0
|
248
|
|
tp@0
|
249 % First, E2E is relative to edge 2
|
|
tp@0
|
250 E2Edist = sqrt( (B5 - zedge1_re2).^2 + redge1_re2.^2 );
|
|
tp@0
|
251
|
|
tp@0
|
252 % B2 will get the coshnyeta values for DF2
|
|
tp@0
|
253 B2 = ( ( B5 - zedge1_re2 ).*( B5 - zR ) + E2Edist.*E2Rdist )./redge1_re2/rR;
|
|
tp@0
|
254 B2 = ( sqrt( B2.^2-1) + B2 ).^nyveclist(2);
|
|
tp@0
|
255 B2 = real( B2 + 1./B2)/2;
|
|
tp@0
|
256
|
|
tp@0
|
257 clear redge1_re2 zedge1_re2
|
|
tp@0
|
258
|
|
tp@0
|
259 %-----------------------------------------------------------------------
|
|
tp@0
|
260 % In directivity function 1, DF1, we need the quantity
|
|
tp@0
|
261 %
|
|
tp@0
|
262 % ( ( ze1_re1 -zs_re1 )*( ze1_re1 - ze2_re1 ) + n*m )/rs/re2_re1
|
|
tp@0
|
263
|
|
tp@0
|
264 % Now, we need E2E relative to edge 1
|
|
tp@0
|
265 E2Edist = sqrt( (B4 - zedge2_re1).^2 + redge2_re1.^2 );
|
|
tp@0
|
266
|
|
tp@0
|
267 % B1 will get the coshnyeta values for DF1
|
|
tp@0
|
268
|
|
tp@0
|
269 B1 = (( B4 - zS ).*( B4 - zedge2_re1 ) + S2Edist.*E2Edist)/rS./redge2_re1;
|
|
tp@0
|
270 B1 = ( real(sqrt( B1.^2-1)) + B1 ).^nyveclist(1);
|
|
tp@0
|
271 B1 = ( B1 + 1./B1)/2;
|
|
tp@0
|
272
|
|
tp@0
|
273 clear zedge2_re1 redge2_re1
|
|
tp@0
|
274
|
|
tp@0
|
275 else
|
|
tp@0
|
276 B1 = 1; B2 = 1;
|
|
tp@0
|
277 end
|
|
tp@0
|
278
|
|
tp@0
|
279 % Calculate DF1.*DF2
|
|
tp@0
|
280
|
|
tp@0
|
281 if pathalongplane == 0
|
|
tp@0
|
282
|
|
tp@0
|
283 % Note that here we use E2E re. edge 1!
|
|
tp@0
|
284 B3 = multfac*nyveclist(1)*nyveclist(2)/16/pi^2*dzvec(1)*dzvec(2)*...
|
|
tp@0
|
285 ( sin(nyveclist(1)*(pi + thetaS + thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi + thetaS + thetaedge2_re1 ))) + ...
|
|
tp@0
|
286 sin(nyveclist(1)*(pi + thetaS - thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi + thetaS - thetaedge2_re1 ))) + ...
|
|
tp@0
|
287 sin(nyveclist(1)*(pi - thetaS + thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi - thetaS + thetaedge2_re1 ))) + ...
|
|
tp@0
|
288 sin(nyveclist(1)*(pi - thetaS - thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi - thetaS - thetaedge2_re1 ))))./S2Edist./E2Edist;
|
|
tp@0
|
289
|
|
tp@0
|
290 B3 = B3.*( sin(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))) + ...
|
|
tp@0
|
291 sin(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))) + ...
|
|
tp@0
|
292 sin(nyveclist(2)*(pi - thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi - thetaedge1_re2 + thetaR))) + ...
|
|
tp@0
|
293 sin(nyveclist(2)*(pi - thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi - thetaedge1_re2 - thetaR))))./E2Rdist;
|
|
tp@0
|
294
|
|
tp@0
|
295 else
|
|
tp@0
|
296
|
|
tp@0
|
297 if thetaS == 0
|
|
tp@0
|
298 B3 = multfac*nyveclist(1)*nyveclist(2)/2/pi^2*dzvec(1)*dzvec(2)*...
|
|
tp@0
|
299 ( sin(nyveclist(1)*(pi + thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi + thetaedge2_re1 ))))./S2Edist./E2Edist;
|
|
tp@0
|
300 else
|
|
tp@0
|
301
|
|
tp@0
|
302 B3 = multfac*nyveclist(1)*nyveclist(2)/4/pi^2*dzvec(1)*dzvec(2)*...
|
|
tp@0
|
303 ( sin(nyveclist(1)*(pi + thetaS + thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi + thetaS + thetaedge2_re1 ))) + ...
|
|
tp@0
|
304 sin(nyveclist(1)*(pi - thetaS + thetaedge2_re1 ))./(B1- cos(nyveclist(1)*(pi - thetaS + thetaedge2_re1 ))))./S2Edist./E2Edist;
|
|
tp@0
|
305 end
|
|
tp@0
|
306
|
|
tp@0
|
307
|
|
tp@0
|
308 B3 = B3.*( sin(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))) + ...
|
|
tp@0
|
309 sin(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))))./E2Rdist;
|
|
tp@0
|
310 end
|
|
tp@0
|
311
|
|
tp@0
|
312 %-------------------------------------------------------------------------
|
|
tp@0
|
313 % Determine in which sample slots the amplitude contributions should
|
|
tp@0
|
314 % be added, based on the total distance.
|
|
tp@0
|
315
|
|
tp@0
|
316 % B2 is number of sample slots as non-integers
|
|
tp@0
|
317
|
|
tp@0
|
318 B2 = ( S2Edist + E2Edist + E2Rdist - R_irstart)/(cair/fs)+1;
|
|
tp@0
|
319 clear S2Edist E2Edist E2Rdist
|
|
tp@0
|
320
|
|
tp@0
|
321 % B1 is the sampleslot numbers in integer number
|
|
tp@0
|
322 B1 = floor(B2);
|
|
tp@0
|
323
|
|
tp@0
|
324 % B2 is a value between 0 and 1 stating how much of dh that should be added
|
|
tp@0
|
325 % to the first sample slot, i.e. the one given on B1. The rest of dh should
|
|
tp@0
|
326 % be added to the following sample slot.
|
|
tp@0
|
327 B2 = B2 - B1;
|
|
tp@0
|
328
|
|
tp@0
|
329 % New approach: simply state that the amplitudevalue, mult. by 1-B2, should
|
|
tp@0
|
330 % be placed in the integer slot given in B1. In addition, the
|
|
tp@0
|
331 % amplitudevalues, mult by B2, should be placed in the integer slot given
|
|
tp@0
|
332 % by B1+1.
|
|
tp@0
|
333
|
|
tp@0
|
334 B1 = [B1;B1+1];
|
|
tp@0
|
335 B3 = [B3.*(1-B2);B3.*B2];
|
|
tp@0
|
336
|
|
tp@0
|
337 ir = zeros(max(B1),1 );
|
|
tp@0
|
338
|
|
tp@0
|
339 [B1,sortvec] = sort(B1);
|
|
tp@0
|
340 B3 = cumsum(B3(sortvec));
|
|
tp@0
|
341 ivstep = (find(diff([(B1);(B1(end))+1])));
|
|
tp@0
|
342 ir(B1(ivstep)) = diff([0;B3(ivstep)]);
|
|
tp@0
|
343
|
|
tp@0
|
344 % Here we could have the possibility to save space by cutting out initial zeros and
|
|
tp@0
|
345 % set the value ninit correspondingly to the number of removed zeros. However, the rest of
|
|
tp@0
|
346 % the EDB1 functions have not implented support for this.
|
|
tp@0
|
347
|
|
tp@0
|
348 ninit = 0;
|
