|
tp@0
|
1 function [ir,initdelay,singularterm] = EDB1wedge1st_int(fs,closwedang,rs,thetas,zs,rr,thetar,zr,zw,Method,R_irstart,bc)
|
|
tp@0
|
2 % EDB1wedge1st_int - Gives the 1st order diffraction impulse response.
|
|
tp@0
|
3 % Gives the 1st order diffraction impulse response
|
|
tp@0
|
4 % for a point source irradiating a finite wedge. A free-field impulse response amplitude
|
|
tp@0
|
5 % of 1/r is used. An accurate integration method is used
|
|
tp@0
|
6 % so receivers can be placed right at the zone boundaries.
|
|
tp@0
|
7 %
|
|
tp@0
|
8 % Output parameters:
|
|
tp@0
|
9 % ir the impulse response
|
|
tp@0
|
10 % initdelay the excluded initial time delay, from the point source switch-on until the start of ir.
|
|
tp@0
|
11 % singularterm a vector, [1 4], of 0 and 1 that indicates if one or more of the four
|
|
tp@0
|
12 % beta terms is singular, i.e., right on a zone boundary.
|
|
tp@0
|
13 % If there is one or more elements with the value 1, the
|
|
tp@0
|
14 % direct sound or a specular reflection needs to be given
|
|
tp@0
|
15 % half its specular value.
|
|
tp@0
|
16 %
|
|
tp@0
|
17 % Input parameters:
|
|
tp@0
|
18 % fs sampling frequency
|
|
tp@0
|
19 % closwedang closed wedge angle
|
|
tp@0
|
20 % rs, thetas, zs cyl. coordinates of the source pos. (zs = 0)
|
|
tp@0
|
21 % rr, thetar, zr cyl. coordinates of the receiver pos
|
|
tp@0
|
22 % zw the edge ends given as [z1 z2]
|
|
tp@0
|
23 % Method 'New' (the new method by Svensson-Andersson-Vanderkooy) or 'Van'
|
|
tp@0
|
24 % (Vanderkooys old method) or 'KDA' (Kirchhoff approximation).
|
|
tp@0
|
25 % R_irstart (optional) if this is included, the impulse response time zero will
|
|
tp@0
|
26 % correspond to this distance. Otherwise, the time zero
|
|
tp@0
|
27 % corresponds to the distance R0 via the apex point.
|
|
tp@0
|
28 % bc (optional) boundary condition, [1 1] (default) => rigid-rigid, [-1 -1] => soft-soft
|
|
tp@0
|
29 % [1 -1] => rigid-soft. First plane should be the ref. plane.
|
|
tp@0
|
30 % CAIR (global) The speed of sound
|
|
tp@0
|
31 %
|
|
tp@0
|
32 % Uses the functions EDB1quadstep and EDB1betaoverml for numerical integration.
|
|
tp@0
|
33 %
|
|
tp@0
|
34 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
35 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
36 %
|
|
tp@0
|
37 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
38 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
39 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
40 %
|
|
tp@0
|
41 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
42 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
43 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
44 %
|
|
tp@0
|
45 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
46 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
47 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
48 % Peter Svensson (svensson@iet.ntnu.no) 20060607
|
|
tp@0
|
49 %
|
|
tp@0
|
50 % [ir,initdelay,singularterm] = EDB1wedge1st_int(fs,closwedang,rs,thetas,zs,rr,thetar,zr,zw,Method,R_irstart,bc);
|
|
tp@0
|
51
|
|
tp@0
|
52 global CAIR
|
|
tp@0
|
53 localshowtext = 0;
|
|
tp@0
|
54
|
|
tp@0
|
55 if localshowtext
|
|
tp@0
|
56 disp(' ')
|
|
tp@0
|
57 disp('Reference solution for wedge 1st')
|
|
tp@0
|
58 end
|
|
tp@0
|
59
|
|
tp@0
|
60 if nargin < 11
|
|
tp@0
|
61 R_irstart = 0;
|
|
tp@0
|
62 end
|
|
tp@0
|
63 if nargin < 12
|
|
tp@0
|
64 bc = [1 1];
|
|
tp@0
|
65 else
|
|
tp@0
|
66 if bc(1)*bc(2)~=1 & bc(1) ~= 1
|
|
tp@0
|
67 error('ERROR: Only all-rigid wedges are implemented');
|
|
tp@0
|
68 end
|
|
tp@0
|
69 end
|
|
tp@0
|
70
|
|
tp@0
|
71 %------------------------------------------------------------------
|
|
tp@0
|
72 % If the wedge has the length zero, return a zero IR immediately
|
|
tp@0
|
73
|
|
tp@0
|
74 if zw(2)-zw(1) == 0
|
|
tp@0
|
75 ir = [0 0].';
|
|
tp@0
|
76 initdelay = 0;
|
|
tp@0
|
77 singularterm = [0 0 0 0];
|
|
tp@0
|
78 return
|
|
tp@0
|
79 end
|
|
tp@0
|
80
|
|
tp@0
|
81 %----------------------------------------------
|
|
tp@0
|
82 % Method + Some edge constants
|
|
tp@0
|
83
|
|
tp@0
|
84 ny = pi/(2*pi-closwedang);
|
|
tp@0
|
85
|
|
tp@0
|
86 initdelay = R_irstart/CAIR;
|
|
tp@0
|
87 sampconst = CAIR/fs;
|
|
tp@0
|
88
|
|
tp@0
|
89 if Method(1) == 'N' | Method(1) == 'n'
|
|
tp@0
|
90 Method = 'New';
|
|
tp@0
|
91 else
|
|
tp@0
|
92 error(['ERROR: The method ',Method,' has not been implemented yet!']);
|
|
tp@0
|
93 end
|
|
tp@0
|
94
|
|
tp@0
|
95 tol = 1e-11; % The relative accuracy for the numerical integration
|
|
tp@0
|
96 % It was found that 1e-6 generated a
|
|
tp@0
|
97 % substantial error for some cases (PC
|
|
tp@0
|
98 % 041129)
|
|
tp@0
|
99
|
|
tp@0
|
100 zrelmax = 0.1; % The part of the edge that should use the analyt. approx.
|
|
tp@0
|
101
|
|
tp@0
|
102 za = (rs*zr+rr*zs)/(rs+rr); % The apex point: the point on the edge with the shortest path
|
|
tp@0
|
103
|
|
tp@0
|
104 % Move the z-values so that za ends up in z = 0.
|
|
tp@0
|
105
|
|
tp@0
|
106 Dz = abs(zs-zr);
|
|
tp@0
|
107 zs = zs - za;
|
|
tp@0
|
108 zr = zr - za;
|
|
tp@0
|
109 zw = zw - za; % zw is a two-element vector
|
|
tp@0
|
110 za = 0;
|
|
tp@0
|
111
|
|
tp@0
|
112 rs2 = rs^2;
|
|
tp@0
|
113 rr2 = rr^2;
|
|
tp@0
|
114
|
|
tp@0
|
115 R0 = sqrt( (rs+rr)^2 + (zr-zs)^2 );
|
|
tp@0
|
116
|
|
tp@0
|
117 % Calculate the sinnyfi and cosnyfi values for use in the four beta-terms
|
|
tp@0
|
118 % later
|
|
tp@0
|
119
|
|
tp@0
|
120 fivec = pi + thetas*[1 1 -1 -1] + thetar*[1 -1 1 -1];
|
|
tp@0
|
121
|
|
tp@0
|
122 absnyfivec = abs(ny*fivec);
|
|
tp@0
|
123
|
|
tp@0
|
124 sinnyfivec = sin(ny*fivec);
|
|
tp@0
|
125 cosnyfivec = cos(ny*fivec);
|
|
tp@0
|
126
|
|
tp@0
|
127 if localshowtext
|
|
tp@0
|
128 disp(' ')
|
|
tp@0
|
129 disp([' absnyfivec = ',num2str(absnyfivec)])
|
|
tp@0
|
130 disp([' cosnyfivec = ',num2str(cosnyfivec)])
|
|
tp@0
|
131 end
|
|
tp@0
|
132
|
|
tp@0
|
133 %------------------------------------------------------------------
|
|
tp@0
|
134 % Check which category we have.
|
|
tp@0
|
135 %
|
|
tp@0
|
136 % 1. Is the apex point inside the physical wedge or not? We want to use an
|
|
tp@0
|
137 % analytic approximation for a small part of the apex section (and
|
|
tp@0
|
138 % possibly an ordinary numerical integration for the rest of the apex
|
|
tp@0
|
139 % section. The ordinary numerical integration will be needed if the apex
|
|
tp@0
|
140 % section extends further than what allows the analytic approximation to
|
|
tp@0
|
141 % be used.)
|
|
tp@0
|
142 %
|
|
tp@0
|
143 % apexincluded 0 or 1
|
|
tp@0
|
144 %
|
|
tp@0
|
145 % 2. Will the IR have a step or not? A step is caused if the path lengths
|
|
tp@0
|
146 % to the two edge end points are different. If these two path lengths
|
|
tp@0
|
147 % are different, then we say that the IR has a "tail". This tail has
|
|
tp@0
|
148 % the amplitude 0.5 times the infinite wedge-IR, because there is no
|
|
tp@0
|
149 % no corresponding wedge part on the opposite side of the apex point.
|
|
tp@0
|
150 % We will use an ordinary numerical integration for the tail section.
|
|
tp@0
|
151 %
|
|
tp@0
|
152 % tailincluded 0 or 1
|
|
tp@0
|
153 %
|
|
tp@0
|
154 % Note that if at least one of apexincluded and tailincluded must be 1.
|
|
tp@0
|
155 %
|
|
tp@0
|
156 % In addition to these two choices, we also want to check if we have a
|
|
tp@0
|
157 % perpendicular case or a symmetrical case because these two cases can use
|
|
tp@0
|
158 % a slightly simpler anaytical approximation.
|
|
tp@0
|
159 %
|
|
tp@0
|
160 % We call it a perpendicular case when zs = zr (which can happen only if
|
|
tp@0
|
161 % apexincluded = 1!)
|
|
tp@0
|
162 %
|
|
tp@0
|
163 % perpcase 0 or 1
|
|
tp@0
|
164 %
|
|
tp@0
|
165 % We call it a symmetrical case when rs = rr (which could also be a
|
|
tp@0
|
166 % perpendicular case. A symmetrical case might not necessarily
|
|
tp@0
|
167 % include the apex point. If the apex point is not included then we are not
|
|
tp@0
|
168 % really interested in whether we have a symmetrical case or not since the
|
|
tp@0
|
169 % analytic approximation is used only for the apex section!)
|
|
tp@0
|
170 %
|
|
tp@0
|
171 % symmcase 0 or 1
|
|
tp@0
|
172 %
|
|
tp@0
|
173 % Since both the perpendicular case and the symmetric case can use the same
|
|
tp@0
|
174 % simpler analytic approximation, we also denote those cases that are
|
|
tp@0
|
175 % neither "skew cases". Again, this is relevant only when the apex is
|
|
tp@0
|
176 % included.
|
|
tp@0
|
177 %
|
|
tp@0
|
178 % skewcase 0 or 1
|
|
tp@0
|
179
|
|
tp@0
|
180 perpcase = 0;
|
|
tp@0
|
181 if zs == zr
|
|
tp@0
|
182 perpcase = 1;
|
|
tp@0
|
183 end
|
|
tp@0
|
184
|
|
tp@0
|
185 symmcase = 0;
|
|
tp@0
|
186 if rs == rr
|
|
tp@0
|
187 symmcase = 1;
|
|
tp@0
|
188 end
|
|
tp@0
|
189
|
|
tp@0
|
190 skewcase = 0;
|
|
tp@0
|
191 if perpcase == 0 & symmcase == 0
|
|
tp@0
|
192 skewcase = 1;
|
|
tp@0
|
193 end
|
|
tp@0
|
194
|
|
tp@0
|
195 apexincluded = 1;
|
|
tp@0
|
196 if sign( zw(1)*zw(2) ) == 1
|
|
tp@0
|
197 apexincluded = 0;
|
|
tp@0
|
198 end
|
|
tp@0
|
199
|
|
tp@0
|
200 Rendnear = sqrt( rs2 + (zw(1)-zs)^2 ) + sqrt( rr2 + (zw(1)-zr)^2 );
|
|
tp@0
|
201 Rendfar = sqrt( rs2 + (zw(2)-zs)^2 ) + sqrt( rr2 + (zw(2)-zr)^2 );
|
|
tp@0
|
202
|
|
tp@0
|
203 if Rendnear > Rendfar
|
|
tp@0
|
204 temp = Rendnear;
|
|
tp@0
|
205 Rendnear = Rendfar;
|
|
tp@0
|
206 Rendfar = temp;
|
|
tp@0
|
207 end
|
|
tp@0
|
208
|
|
tp@0
|
209 tailincluded = 1;
|
|
tp@0
|
210 if abs(Rendfar-Rendnear) < eps*1e6
|
|
tp@0
|
211 tailincluded = 0;
|
|
tp@0
|
212 end
|
|
tp@0
|
213
|
|
tp@0
|
214 %------------------------------------------------------------------
|
|
tp@0
|
215 % Find out where the IR should start and end
|
|
tp@0
|
216
|
|
tp@0
|
217 arrivalsampnumb = round((R0-R_irstart)/sampconst);
|
|
tp@0
|
218 endsampnumb_apex = round((Rendnear-R_irstart)/sampconst);
|
|
tp@0
|
219 if tailincluded == 1
|
|
tp@0
|
220 endsampnumb_tail = round((Rendfar-R_irstart)/sampconst);
|
|
tp@0
|
221 ir = zeros(endsampnumb_tail+1,1);
|
|
tp@0
|
222 else
|
|
tp@0
|
223 ir = zeros(endsampnumb_apex+1,1);
|
|
tp@0
|
224 end
|
|
tp@0
|
225
|
|
tp@0
|
226 %------------------------------------------------------------------
|
|
tp@0
|
227 % Construct vectors with the integration intervals along the edge
|
|
tp@0
|
228 %
|
|
tp@0
|
229 % For the apex section, which will include one part along the negative
|
|
tp@0
|
230 % branch and one part along the positive branch of the wedge
|
|
tp@0
|
231 % we carry out the integration only along the positive branch and multiply
|
|
tp@0
|
232 % by two.
|
|
tp@0
|
233 % For the tail section of the wedge (if it is included), we will also carry
|
|
tp@0
|
234 % out the numerical intregration along the positive branch - even if the
|
|
tp@0
|
235 % physical wedge actually has its tail secction along the negative branch.
|
|
tp@0
|
236 %
|
|
tp@0
|
237 % This choice implies means that we must always find the positive solutions
|
|
tp@0
|
238 % to the equation where we find z-values along the wedge that correspond to
|
|
tp@0
|
239 % the integration time values.
|
|
tp@0
|
240
|
|
tp@0
|
241 % Calculate the integration intervals for the apex section
|
|
tp@0
|
242
|
|
tp@0
|
243 if apexincluded == 1
|
|
tp@0
|
244
|
|
tp@0
|
245 % x is the path length (m+l) which corresponds to the end of each
|
|
tp@0
|
246 % integration interval = time sample.
|
|
tp@0
|
247 % This is the part of the edge which is symmetrical around the apex point.
|
|
tp@0
|
248
|
|
tp@0
|
249 x = sampconst*(arrivalsampnumb+[0.5:1:(endsampnumb_apex-arrivalsampnumb)+0.5]) + R_irstart;
|
|
tp@0
|
250 nx = length(x);
|
|
tp@0
|
251
|
|
tp@0
|
252 if x(nx) > Rendnear
|
|
tp@0
|
253 x(nx) = Rendnear;
|
|
tp@0
|
254 if nx > 1
|
|
tp@0
|
255 if x(nx) < x(nx-1)
|
|
tp@0
|
256 error('ERROR: Strange error number 1')
|
|
tp@0
|
257 end
|
|
tp@0
|
258 end
|
|
tp@0
|
259 end
|
|
tp@0
|
260 x2 = x.^2;
|
|
tp@0
|
261 % ptc note: x is the same as L which is used in other derivations
|
|
tp@0
|
262
|
|
tp@0
|
263 % Convert the x-values to z-values along the edge. These z-values will
|
|
tp@0
|
264 % define the integration intervals.
|
|
tp@0
|
265
|
|
tp@0
|
266 % ptc new code
|
|
tp@0
|
267 M02 = rs*rs + zs*zs;
|
|
tp@0
|
268 L02 = rr*rr + zr*zr;
|
|
tp@0
|
269 K = M02 - L02 - x2;
|
|
tp@0
|
270
|
|
tp@0
|
271 diffz = (zs - zr);
|
|
tp@0
|
272 denom = diffz^2 - x2;
|
|
tp@0
|
273 a = (2*x2*zr - K*diffz) ./ denom;
|
|
tp@0
|
274 b = ((K.^2 / 4) - x2*L02) ./ denom;
|
|
tp@0
|
275
|
|
tp@0
|
276 zrange_apex = -a/2 + sqrt((a.^2)/4 - b);
|
|
tp@0
|
277
|
|
tp@0
|
278 % Finally the endpoint of the integral of the analytical approximation
|
|
tp@0
|
279 % (=zrangespecial) should be only up to z = 0.01
|
|
tp@0
|
280 % or whatever is specified as zrelmax for the analytic approximation
|
|
tp@0
|
281
|
|
tp@0
|
282 zrangespecial = zrelmax*min([rs,rr]);
|
|
tp@0
|
283 % The explicit values below were just tested to get identical results with
|
|
tp@0
|
284 % another implementation
|
|
tp@0
|
285 % zrangespecial = 0.01629082265142;
|
|
tp@0
|
286 % zrangespecial = 0.01530922080939; %0.00521457034952;
|
|
tp@0
|
287 % zrangespecial = 0.00510940311841;
|
|
tp@0
|
288 % disp('End of apex range = :')
|
|
tp@0
|
289 % zrange_apex(1)
|
|
tp@0
|
290
|
|
tp@0
|
291 if zrangespecial > zrange_apex(1)
|
|
tp@0
|
292 zrangespecial = abs(zrange_apex(1));
|
|
tp@0
|
293 splitinteg = 0;
|
|
tp@0
|
294 else
|
|
tp@0
|
295 splitinteg = 1;
|
|
tp@0
|
296 end
|
|
tp@0
|
297 end
|
|
tp@0
|
298
|
|
tp@0
|
299 % Calculate the integration intervals for the tail section
|
|
tp@0
|
300
|
|
tp@0
|
301 if tailincluded == 1
|
|
tp@0
|
302
|
|
tp@0
|
303 % x is the path length (m+l) which corresponds to the end of each
|
|
tp@0
|
304 % integration interval = time sample.
|
|
tp@0
|
305 % This is the part of the edge which is outside the symmetrical part around the apex point.
|
|
tp@0
|
306
|
|
tp@0
|
307 x = sampconst*(endsampnumb_apex+[-0.5:1:(endsampnumb_tail-endsampnumb_apex)+0.5]) + R_irstart;
|
|
tp@0
|
308 nx = length(x);
|
|
tp@0
|
309 if x(nx) > Rendfar
|
|
tp@0
|
310 x(nx) = Rendfar;
|
|
tp@0
|
311 if nx > 1
|
|
tp@0
|
312 if x(nx) < x(nx-1)
|
|
tp@0
|
313 error('ERROR: Strange error number 2')
|
|
tp@0
|
314 end
|
|
tp@0
|
315 end
|
|
tp@0
|
316 end
|
|
tp@0
|
317 if x(1) < Rendnear
|
|
tp@0
|
318 x(1) = Rendnear;
|
|
tp@0
|
319 if nx > 1
|
|
tp@0
|
320 if x(2) < x(1)
|
|
tp@0
|
321 error('ERROR: Strange error number 2')
|
|
tp@0
|
322 end
|
|
tp@0
|
323 end
|
|
tp@0
|
324 end
|
|
tp@0
|
325 x2 = x.^2;
|
|
tp@0
|
326
|
|
tp@0
|
327 % ptc new code
|
|
tp@0
|
328 M02 = rs*rs + zs*zs;
|
|
tp@0
|
329 L02 = rr*rr + zr*zr;
|
|
tp@0
|
330 K = M02 - L02 - x2;
|
|
tp@0
|
331 diffz = (zs - zr);
|
|
tp@0
|
332 denom = diffz^2 - x2;
|
|
tp@0
|
333 a = (2*x2*zr - K*diffz) ./ denom;
|
|
tp@0
|
334 b = ((K.^2 / 4) - x2*L02) ./ denom;
|
|
tp@0
|
335
|
|
tp@0
|
336 zrange_tail = -a/2 + sqrt((a.^2)/4 - b);
|
|
tp@0
|
337 end
|
|
tp@0
|
338
|
|
tp@0
|
339 if localshowtext
|
|
tp@0
|
340 if apexincluded
|
|
tp@0
|
341 if tailincluded
|
|
tp@0
|
342 disp([int2str(length(zrange_apex)+length(zrange_tail)),' elements'])
|
|
tp@0
|
343 else
|
|
tp@0
|
344 disp([int2str(length(zrange_apex)),' elements'])
|
|
tp@0
|
345 end
|
|
tp@0
|
346 else
|
|
tp@0
|
347 disp([int2str(length(zrange_tail)),' elements'])
|
|
tp@0
|
348 end
|
|
tp@0
|
349 end
|
|
tp@0
|
350
|
|
tp@0
|
351 %----------------------------------------------------------------
|
|
tp@0
|
352 % Compute the impulse responses by carrying out a numerical integration for
|
|
tp@0
|
353 % the little segments that represent each sample slot.
|
|
tp@0
|
354 %
|
|
tp@0
|
355 % For the apex section, there is one little part of the first sample slot
|
|
tp@0
|
356 % which can employ an explicit integration based on the analytical
|
|
tp@0
|
357 % approximation. Often, part of the first sample slot is done with the
|
|
tp@0
|
358 % explicit integration value and another part is done with usual
|
|
tp@0
|
359 % integration.
|
|
tp@0
|
360
|
|
tp@0
|
361 singularterm = [0 0 0 0];
|
|
tp@0
|
362
|
|
tp@0
|
363 if apexincluded == 1
|
|
tp@0
|
364
|
|
tp@0
|
365 %-----------------------------------------------------------
|
|
tp@0
|
366 % First the analytical approximation
|
|
tp@0
|
367
|
|
tp@0
|
368 rho = rr/rs;
|
|
tp@0
|
369 sinpsi = (rs+rr)/R0;
|
|
tp@0
|
370 cospsi = (zr-zs)/R0;
|
|
tp@0
|
371 if localshowtext
|
|
tp@0
|
372 disp([' cospsi = ',num2str(cospsi)])
|
|
tp@0
|
373 end
|
|
tp@0
|
374
|
|
tp@0
|
375 tempfact = (1+rho)^2*sinpsi^2 - 2*rho;
|
|
tp@0
|
376
|
|
tp@0
|
377 useserialexp1 = absnyfivec < 0.01;
|
|
tp@0
|
378 useserialexp2 = abs(absnyfivec - 2*pi) < 0.01;
|
|
tp@0
|
379 useserialexp = useserialexp1 | useserialexp2;
|
|
tp@0
|
380
|
|
tp@0
|
381 singularterm = absnyfivec < 10*eps | abs(absnyfivec - 2*pi) < 10*eps;
|
|
tp@0
|
382 if any(singularterm) & localshowtext
|
|
tp@0
|
383 disp([' Singularity for term ',int2str(find(singularterm))])
|
|
tp@0
|
384 end
|
|
tp@0
|
385
|
|
tp@0
|
386 sqrtB1vec = ( sqrt( 2*(1-cosnyfivec) )*R0*rho/(1+rho)^2/ny ).*(useserialexp==0) + ...
|
|
tp@0
|
387 ( fivec.*sqrt(1-(ny*fivec).^2/12)*R0*rho/(1+rho)^2 ).*(useserialexp1==1) + ...
|
|
tp@0
|
388 ( (fivec-2*pi/ny).*sqrt(1-(ny*fivec-2*pi).^2/12)*R0*rho/(1+rho)^2 ).*(useserialexp2==1);
|
|
tp@0
|
389
|
|
tp@0
|
390 if skewcase == 0, % Use the slightly simpler formulation of the analytical approx.
|
|
tp@0
|
391
|
|
tp@0
|
392 sqrtB3 = sqrt(2)*R0*rho/(1+rho)/sqrt(tempfact);
|
|
tp@0
|
393
|
|
tp@0
|
394 usespecialcase = abs(sqrtB3 - sqrtB1vec) < 1e-14;
|
|
tp@0
|
395
|
|
tp@0
|
396 fifactvec = ( (1-cosnyfivec)./ny^2 ).*(useserialexp==0) + ...
|
|
tp@0
|
397 ( fivec.^2/2.*(1-(ny*fivec).^2/12) ).*(useserialexp1==1) + ...
|
|
tp@0
|
398 ( (fivec-2*pi/ny).^2/2.*(1-(ny*fivec-2*pi).^2/12) ).*(useserialexp2==1);
|
|
tp@0
|
399
|
|
tp@0
|
400 temp1vec = sinnyfivec./( (1+rho)^2 - tempfact.*fifactvec + eps*10);
|
|
tp@0
|
401
|
|
tp@0
|
402
|
|
tp@0
|
403 temp1_2vec = (sinnyfivec+ 10*eps)./( (1+rho)^2 - tempfact.*fifactvec)./(sqrtB1vec+10*eps).*atan( zrangespecial./(sqrtB1vec+eps) );
|
|
tp@0
|
404
|
|
tp@0
|
405 temp3vec = -1./sqrtB3.*atan( zrangespecial./sqrtB3 );
|
|
tp@0
|
406
|
|
tp@0
|
407 approxintvalvec = 2/ny^2*rho*(temp1_2vec + temp1vec*temp3vec);
|
|
tp@0
|
408 approxintvalvec = approxintvalvec.*(sinnyfivec~=0).*(usespecialcase==0);
|
|
tp@0
|
409
|
|
tp@0
|
410 if any(usespecialcase)
|
|
tp@0
|
411 disp('SPECIAL CASE!')
|
|
tp@0
|
412 specialcasevalue = 1/(2*sqrtB3*sqrtB3).*( zrangespecial./(zrangespecial^2 + sqrtB3*sqrtB3) + 1./sqrtB3.*atan(zrangespecial./(sqrtB3+eps)) );
|
|
tp@0
|
413 approxintvalvec = approxintvalvec + 4*R0*R0*rho*rho*rho*sinnyfivec/ny/ny/(1+rho)^4./tempfact.*specialcasevalue.*(usespecialcase==1);
|
|
tp@0
|
414 end
|
|
tp@0
|
415
|
|
tp@0
|
416
|
|
tp@0
|
417 else % Use the more general formulation of the analytical approx.
|
|
tp@0
|
418 B3 = 2*R0*R0*rho*rho/(1+rho)/(1+rho)/tempfact;
|
|
tp@0
|
419 B1vec = sqrtB1vec.^2;
|
|
tp@0
|
420 B2 = -2*R0*(1-rho)*rho*cospsi/(1+rho)/tempfact;
|
|
tp@0
|
421 E1vec = 4*R0^2*rho^3*sinnyfivec./ny^2/(1+rho)^4./tempfact;
|
|
tp@0
|
422 % Error found 050117 by PS. The line below is wrong and should be
|
|
tp@0
|
423 % replaced by the next one.
|
|
tp@0
|
424 % Wrong version: multfact = E1vec.*B2./(B1vec*B2^2 + B1vec.^2 - B3.^2);
|
|
tp@0
|
425 multfact = -E1vec.*B2./(B1vec*B2^2 + (B1vec - B3).^2);
|
|
tp@0
|
426
|
|
tp@0
|
427 % Error found 050118 by PS: the abs in the P1 equation were missing!
|
|
tp@0
|
428 P1 = 0.5*log( abs(B3)*abs(zrangespecial^2 + B1vec)./abs(B1vec)./abs(zrangespecial^2 + B2*zrangespecial + B3 ) );
|
|
tp@0
|
429 P2 = (B1vec-B3)./(sqrtB1vec+10*eps)./B2.*atan(zrangespecial./(sqrtB1vec+10*eps));
|
|
tp@0
|
430 q = 4*B3-B2^2;
|
|
tp@0
|
431
|
|
tp@0
|
432 if q > 0
|
|
tp@0
|
433 sqrtq = sqrt(q);
|
|
tp@0
|
434 F = 2/sqrtq*( atan( (2*zrangespecial+B2)/sqrtq ) - atan( B2/sqrtq ) );
|
|
tp@0
|
435 elseif q < 0
|
|
tp@0
|
436 sqrtminq = sqrt(-q);
|
|
tp@0
|
437 % Error found 050118 by PS: the abs in the F equation were missing!
|
|
tp@0
|
438 F = 1./sqrtminq.*log( abs(2*zrangespecial+B2-sqrtminq).*abs(B2+sqrtminq)./abs(2*zrangespecial+B2+sqrtminq)./abs(B2-sqrtminq) );
|
|
tp@0
|
439 else % q = 0
|
|
tp@0
|
440 F = 4*zrangespecial./B2./(2*zrangespecial+B2);
|
|
tp@0
|
441 end
|
|
tp@0
|
442 P3 = (2*(B3-B1vec)-B2^2)/2/B2.*F;
|
|
tp@0
|
443
|
|
tp@0
|
444 % Error found 050118 by PC. The line below is wrong and should be replaced
|
|
tp@0
|
445 % by the one after.
|
|
tp@0
|
446 % Wrong version: approxintvalvec = sum(multfact.*(P1+P2+P3))
|
|
tp@0
|
447 approxintvalvec = multfact.*(P1+P2+P3);
|
|
tp@0
|
448 end
|
|
tp@0
|
449
|
|
tp@0
|
450 approxintvalvec = sum(approxintvalvec.*(1-singularterm));
|
|
tp@0
|
451 if localshowtext
|
|
tp@0
|
452 disp([' Analyt. int = (final IR value) ',num2str(sum(approxintvalvec*(-ny/2/pi)),15)])
|
|
tp@0
|
453 end
|
|
tp@0
|
454
|
|
tp@0
|
455 %-----------------------------------------------------------
|
|
tp@0
|
456 % Numerical integration for the rest of the apex section
|
|
tp@0
|
457
|
|
tp@0
|
458 % The first IR sample will either have part analytic, part numerical
|
|
tp@0
|
459 % integration, or just analytic
|
|
tp@0
|
460
|
|
tp@0
|
461 ir(1+arrivalsampnumb) = ir(1+arrivalsampnumb) + approxintvalvec;
|
|
tp@0
|
462
|
|
tp@0
|
463 if splitinteg == 1
|
|
tp@0
|
464 h = 0.13579*(zrange_apex(1)-zrangespecial);
|
|
tp@0
|
465 x = [zrangespecial zrangespecial+h zrangespecial+2*h (zrangespecial+zrange_apex(1))/2 zrange_apex(1)-2*h zrange_apex(1)-h zrange_apex(1)];
|
|
tp@0
|
466 y = EDB1betaoverml(x,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
467
|
|
tp@0
|
468 Q = [0 0 0];
|
|
tp@0
|
469 Q(:,1) = EDB1quadstep(x(:,1),x(:,3),y(:,1),y(:,2),y(:,3),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
470 Q(:,2) = EDB1quadstep(x(:,3),x(:,5),y(:,3),y(:,4),y(:,5),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
471 Q(:,3) = EDB1quadstep(x(:,5),x(:,7),y(:,5),y(:,6),y(:,7),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
472 ir(1+arrivalsampnumb) = ir(1+arrivalsampnumb) + sum(Q);
|
|
tp@0
|
473
|
|
tp@0
|
474 if localshowtext
|
|
tp@0
|
475 disp([' Split: num value = (final IR value) ',num2str(sum(Q)*(-ny/2/pi)),15])
|
|
tp@0
|
476 disp([' Arrival sample = ',int2str(1+arrivalsampnumb)])
|
|
tp@0
|
477 end
|
|
tp@0
|
478 end
|
|
tp@0
|
479 if localshowtext
|
|
tp@0
|
480 disp([' First IR value = ',num2str(ir(1+arrivalsampnumb)*(-ny/2/pi),15)])
|
|
tp@0
|
481 disp(' ')
|
|
tp@0
|
482 end
|
|
tp@0
|
483
|
|
tp@0
|
484 % The numerical integration for the rest of the apex section
|
|
tp@0
|
485
|
|
tp@0
|
486 zstarts = zrange_apex(1:end-1).';
|
|
tp@0
|
487 zends = zrange_apex(2:end).';
|
|
tp@0
|
488
|
|
tp@0
|
489 h = 0.13579*(zends-zstarts);
|
|
tp@0
|
490 nslots = length(h);
|
|
tp@0
|
491 x = [zstarts zstarts+h zstarts+2*h (zstarts+zends)/2 zends-2*h zends-h zends];
|
|
tp@0
|
492 x = reshape(x,7*nslots,1);
|
|
tp@0
|
493 y = EDB1betaoverml(x,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
494 y = reshape(y,nslots,7);
|
|
tp@0
|
495 x = reshape(x,nslots,7);
|
|
tp@0
|
496 Q = zeros(nslots,3);
|
|
tp@0
|
497
|
|
tp@0
|
498 Q(:,1) = EDB1quadstep(x(:,1),x(:,3),y(:,1),y(:,2),y(:,3),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
499 Q(:,2) = EDB1quadstep(x(:,3),x(:,5),y(:,3),y(:,4),y(:,5),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
500 Q(:,3) = EDB1quadstep(x(:,5),x(:,7),y(:,5),y(:,6),y(:,7),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
501 ir(arrivalsampnumb+[2:length(zrange_apex)]) = ir(arrivalsampnumb+[2:length(zrange_apex)]) + sum(Q.').';
|
|
tp@0
|
502 ir = ir*(-ny/2/pi); % Mult by 2 because analyt int. is half the wedge
|
|
tp@0
|
503
|
|
tp@0
|
504 end
|
|
tp@0
|
505
|
|
tp@0
|
506 %--------------------------------------------------------------------------
|
|
tp@0
|
507 % If we have a non-symmetrical case, then we must integrate for the long
|
|
tp@0
|
508 % end of the edge as well.
|
|
tp@0
|
509
|
|
tp@0
|
510 if tailincluded == 1
|
|
tp@0
|
511 zstarts = zrange_tail(1:end-1).';
|
|
tp@0
|
512 zends = zrange_tail(2:end).';
|
|
tp@0
|
513
|
|
tp@0
|
514 h = 0.13579*(zends-zstarts);
|
|
tp@0
|
515 nslots = length(h);
|
|
tp@0
|
516 x = [zstarts zstarts+h zstarts+2*h (zstarts+zends)/2 zends-2*h zends-h zends];
|
|
tp@0
|
517 x = reshape(x,7*nslots,1);
|
|
tp@0
|
518 y = EDB1betaoverml(x,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
519 y = reshape(y,nslots,7);
|
|
tp@0
|
520 x = reshape(x,nslots,7);
|
|
tp@0
|
521 Q = zeros(nslots,3);
|
|
tp@0
|
522 Q(:,1) = EDB1quadstep(x(:,1),x(:,3),y(:,1),y(:,2),y(:,3),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
523 Q(:,2) = EDB1quadstep(x(:,3),x(:,5),y(:,3),y(:,4),y(:,5),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
524 Q(:,3) = EDB1quadstep(x(:,5),x(:,7),y(:,5),y(:,6),y(:,7),tol,rs,rr,zs,zr,ny,sinnyfivec,cosnyfivec);
|
|
tp@0
|
525
|
|
tp@0
|
526 ir(endsampnumb_apex+[1:length(zrange_tail)-1]) = ir(endsampnumb_apex+[1:length(zrange_tail)-1]) + (sum(Q.').')*(-ny/4/pi);
|
|
tp@0
|
527
|
|
tp@0
|
528 end
|
|
tp@0
|
529
|
|
tp@0
|
530 ir = real(ir);
|
|
tp@0
|
531
|
|
tp@0
|
532 if localshowtext
|
|
tp@0
|
533 disp([' and finally,'])
|
|
tp@0
|
534 disp([' the first IR values = ',num2str(ir(0+arrivalsampnumb),15)])
|
|
tp@0
|
535 disp([' ',num2str(ir(1+arrivalsampnumb),15)])
|
|
tp@0
|
536 disp([' ',num2str(ir(2+arrivalsampnumb),15)])
|
|
tp@0
|
537 disp(' ')
|
|
tp@0
|
538 end
|
