tp@0: function [ir,ninit] = EDB1wedge2nd(cylS,cylR,cylE2_r1,cylE1_r2,...
tp@0: nyveclist,edgelengthlist,dzvec,method,pathalongplane,R_irstart,bc,cair,fs)
tp@0: % EDB1wedge2nd - Gives the 2nd-order diffraction IR.
tp@0: %
tp@0: % NB!!! Only edges that are in-plane with each other can be handled with
tp@0: % this version. A small extension is needed to handle non-in-plane edge
tp@0: % pairs.
tp@0: %
tp@0: % NB!!! Known limitation/error: geometries like reflector clouds with
tp@0: % several surfaces in the same plane will give some errors: obstruction
tp@0: % checks are not done properly for edge-to-edge pairs if they are in the
tp@0: % same plane. Some removal of non-allowed edge-to-edge pairs is done: if a
tp@0: % path from one edge to another edges starts by crossing the first edge's
tp@0: % own plane, then that edge-to-edge combination is shut off. This could be
tp@0: % expressed as if a second edge is "behind" the first edge. It corresponds
tp@0: % somewhat to the check for the image source method if an image source is
tp@0: % behind or in front of its last potential reflection surface -  a
tp@0: % "validity" check. But, the subsequent "obstruction" check is not done
tp@0: % yet.
tp@0: %
tp@0: %   ir     The 2nd order diffraction impulse response for two
tp@0: %          edges that share one plane.
tp@0: %   ninit  The initial delay that has been cut away before the
tp@0: %          beginning of ir. Note that the time zero of the ir (including ninit)
tp@0: %		   corresponds to R_irstart.
tp@0: %
tp@0: % with input parameters
tp@0: %	cylS	vector containing [rS,thetaS,zS] = cyl. coordinates of the source rel. to the first edge
tp@0: %	cylR	vector containing [rR,thetaR,zR] = cyl. coordinates of the receiver rel. to the second edge
tp@0: %	cylE2_r1	matrix containing [rstart_r1,thetastart_r1,zstart_r1;rend_r1,thetaend_r1,zend_r1]
tp@0: %				= cyl. coordinates of the start and end points of the second edge rel. to the first edge
tp@0: %	cylE1_r2	matrix containing [rstart_r2,thetastart_r2,zstart_r2;rend_r2,thetaend_r2,zend_r2]
tp@0: %				= cyl. coordinates of the start and end points of the first edge rel. to the second edge
tp@0: %   nyveclist
tp@0: %          	vector with the two wedge indices of the two edges
tp@0: %   edgelengthlist 
tp@0: %           Matrix, [2 2], with the start and end values of the two edges:
tp@0: %           Row 1 has the start and end values of edge 1
tp@0: %           Row 2 has the start and end values of edge 2
tp@0: %           For a fully visible/active edge, the start value = 0 and the
tp@0: %           end value equals the length (in meters) of the edge.
tp@0: %	elemsize_2nd
tp@0: %		   	the number of elements per wavelength at the sampling frequency.
tp@0: %		   	Recommended value 0.5, and this value is chosen if no other value
tp@0: %		   	has been set. One can experiment with lower values (0.25, 0.1, ...)
tp@0: %		   	for comparisons - the error increases with frequency but if first
tp@0: %		   	order diffraction and specular response is masking the second order
tp@0: %		   	diffraction, low values of elemsize_2nd might be perfectly acceptable.
tp@0: %	method	'n' for the New method or 'v' for Vanderkooys method
tp@0: %   pathalongplane      1 or 0 depending on if the path runs along a plane
tp@0: %           or not
tp@0: %   R_irstart (optional)
tp@0: %          	a distance to which the ir time zero will correspond
tp@0: %	bc (optional)
tp@0: %			a matrix containing [refl11 refl12;refl21 refl22] where the values are +1 or -1
tp@0: %			indicating the reflection coefficients of the four wedge planes involved.
tp@0: %			+1 means rigid and -1 means soft (pressure-release)
tp@0: %			refl11 is edge 1, plane 1 (refplane)
tp@0: %			refl12 is edge 1, plane 2
tp@0: %			refl21 is edge 2, plane 1 (refplane)
tp@0: %			refl22 is edge 2, plane 2
tp@0: %			default is all-rigid
tp@0: %   CAIR (global)
tp@0: %          	the speed of sound
tp@0: %	FSAMP (global)
tp@0: %			the sampling frequency
tp@0: %   BIGEDGESTEPMATRIX (global)
tp@0: %           a matrix, size [N,2], where N = nedge1*nedge2. nedge1 is the number of edge
tp@0: %           elements that edge1 is sub-divided into and nedge2 is the
tp@0: %           number of edge elements that edge2 is sub-divided into. The
tp@0: %           matrix contains values between 0 and 1, representing all
tp@0: %           combinations of relative positions along the two edges for the
tp@0: %           mid-points of these edge elements. They are organized like
tp@0: %           below, if nedge1 = 10 and nedge2 = 25:
tp@0: %           BIGEDGESTEPMATRIX = 
tp@0: %                  [ 0.05   0.02;
tp@0: %                    0.05   0.06;
tp@0: %                    0.05   0.10;
tp@0: %                    ...
tp@0: %                    0.05   0.98;
tp@0: %                    0.15   0.02;
tp@0: %                    ...
tp@0: %                    0.95   0.98];
tp@0: %           These values are multiplied with the lengths of the respective
tp@0: %           edges to give the actual position, in meters along each edge.
tp@0: %
tp@0: % Uses no special functions
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This file is part of the Edge Diffraction Toolbox by Peter Svensson.                       
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify       
tp@0: %   it under the terms of the GNU General Public License as published by the Free Software     
tp@0: %   Foundation, either version 3 of the License, or (at your option) any later version.        
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is distributed in the hope that it will be useful,       
tp@0: %   but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS  
tp@0: %   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.             
tp@0: %                                                                                              
tp@0: %   You should have received a copy of the GNU General Public License along with the           
tp@0: %   Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20080508
tp@0: %
tp@0: % [ir,ninit] = EDB1wedge2ndnew(cylS,cylR,cylE2_r1,cylE1_r2,...
tp@0: % nyveclist,edgelengthlist,dzvec,method,pathalongplane,BigB,R_irstart,bc,cair,fs);
tp@0: 
tp@0: global BIGEDGESTEPMATRIX
tp@0: 
tp@0: multfac = 1/(pathalongplane + 1);
tp@0: 
tp@0: %-----------------------------------------------------------------------
tp@0: % Extract the individual parameters
tp@0: 
tp@0: rS     = cylS(1);
tp@0: thetaS = cylS(2);
tp@0: zS     = cylS(3);
tp@0: 
tp@0: rR     = cylR(1);
tp@0: thetaR = cylR(2);
tp@0: zR     = cylR(3);
tp@0: 
tp@0: rE1_r1     = cylE2_r1(1,1);
tp@0: thetaE1_r1 = cylE2_r1(1,2);
tp@0: zE1_r1     = cylE2_r1(1,3);
tp@0: rE2_r1     = cylE2_r1(2,1);
tp@0: thetaE2_r1 = cylE2_r1(2,2);
tp@0: zE2_r1     = cylE2_r1(2,3);
tp@0: 
tp@0: rE1_r2     = cylE1_r2(1,1);
tp@0: thetaE1_r2 = cylE1_r2(1,2);
tp@0: zE1_r2     = cylE1_r2(1,3);
tp@0: rE2_r2     = cylE1_r2(2,1);
tp@0: thetaE2_r2 = cylE1_r2(2,2);
tp@0: zE2_r2     = cylE1_r2(2,3);
tp@0: 
tp@0: edgesareinplane  = ( abs(thetaE1_r1 - thetaE2_r1) < 1e-6 );
tp@0: 
tp@0: swapbigmatrix = 0;
tp@0: zmid1 = (edgelengthlist(1,1)+edgelengthlist(1,2))/2;
tp@0: zmid2_re1 = ( zE1_r1 + zE2_r1  )/2;
tp@0: midpointdist = sqrt( rE1_r1.^2  + (zmid2_re1-zmid1).^2); 
tp@0: reledgetoedgedist = midpointdist/max(dzvec);
tp@0: 
tp@0: % % % % % % % % % % % if reledgetoedgedist < 12
tp@0: % % % % % % % % % % %     refinefactor = ceil(12/reledgetoedgedist);
tp@0: % % % % % % % % % % % %      disp(['   Refining by factor ',int2str(refinefactor)])
tp@0: % % % % % % % % % % %     newndivvec = round(edgelengthlist(:,2)./dzvec.').'*refinefactor;
tp@0: % % % % % % % % % % %     ivmatrix = EDB1creindexmatrix(newndivvec);
tp@0: % % % % % % % % % % %     [nedgeelcombs,slask] = size(ivmatrix);
tp@0: % % % % % % % % % % %     BIGE = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),:);
tp@0: % % % % % % % % % % %     swapbigmatrix = 1;
tp@0: % % % % % % % % % % %     dzvec = dzvec/refinefactor;
tp@0: % % % % % % % % % % % end
tp@0: 
tp@0: % E2E, first part
tp@0: 
tp@0: if swapbigmatrix == 0
tp@0:     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:     zedge2_re1 = zE1_r1 + BIGEDGESTEPMATRIX(:,2)*( zE2_r1 - zE1_r1 );     % zedge2_re1 are the z-values of edge2, expressed in edge 1
tp@0:     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:     zedge1_re2 = zE1_r2 + BIGEDGESTEPMATRIX(:,1)*( zE2_r2 - zE1_r2 ); 
tp@0: else
tp@0:     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:     zedge2_re1 = zE1_r1 + BIGE(:,2)*( zE2_r1 - zE1_r1 );     % zedge2_re1 are the z-values of edge2, expressed in edge 1
tp@0:     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:     zedge1_re2 = zE1_r2 + BIGE(:,1)*( zE2_r2 - zE1_r2 ); 
tp@0: end
tp@0: 
tp@0: % S2E
tp@0: S2Edist = sqrt( rS^2 + ( B4 - zS ).^2);
tp@0: 
tp@0: % E2R
tp@0: E2Rdist = sqrt( rR^2 + ( B5 - zR ).^2);
tp@0: 
tp@0: % E2E, second part
tp@0: if swapbigmatrix == 0
tp@0:     if edgesareinplane == 0,        
tp@0:         
tp@0: %         error('ERROR: r and theta calculation for non-in-plane edges not implemented yet!')
tp@0: 
tp@0:         % First we need to reconvert the cylindrical coordinates of
tp@0:         % edge2-re1 into cartesian! Then we can use the EDB1coordtrans.
tp@0:         % We define our own cartesian coord syst such that the reference
tp@0:         % edge has its starting point in [0 0 0], and the x-axis is along
tp@0:         % the reference plane.
tp@0:         
tp@0:         xE1_re1 = rE1_r1*cos(thetaE1_r1);
tp@0:         yE1_re1 = rE1_r1*sin(thetaE1_r1);
tp@0:         xE2_re1 = rE2_r1*cos(thetaE2_r1);
tp@0:         yE2_re1 = rE2_r1*sin(thetaE2_r1);
tp@0:         
tp@0:         xvec_re1 = xE1_re1 + (xE2_re1 - xE1_re1)*BIGEDGESTEPMATRIX(:,2);
tp@0:         yvec_re1 = yE1_re1 + (yE2_re1 - yE1_re1)*BIGEDGESTEPMATRIX(:,2);
tp@0:         zvec_re1 = zE1_r1 + (zE2_r1 - zE1_r1)*BIGEDGESTEPMATRIX(:,2);       
tp@0:         
tp@0:         % NB!!! Below we need to check if really the full lengths of both
tp@0:         % edges should be used. If S or R see only part of the respective
tp@0:         % edge then the edge-to-edge contribution should be "windowed" too.
tp@0:         [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:          
tp@0:         xE1_re2 = rE1_r2*cos(thetaE1_r2);
tp@0:         yE1_re2 = rE1_r2*sin(thetaE1_r2);
tp@0:         xE2_re2 = rE2_r2*cos(thetaE2_r2);
tp@0:         yE2_re2 = rE2_r2*sin(thetaE2_r2);
tp@0:         
tp@0:         xvec_re2 = xE1_re2 + (xE2_re2 - xE1_re2)*BIGEDGESTEPMATRIX(:,1);
tp@0:         yvec_re2 = yE1_re2 + (yE2_re2 - yE1_re2)*BIGEDGESTEPMATRIX(:,1);
tp@0:         zvec_re2 = zE1_r2 + (zE2_r2 - zE1_r2)*BIGEDGESTEPMATRIX(:,1);
tp@0:                 
tp@0:         % NB!!! Below we need to check if really the full lengths of both
tp@0:         % edges should be used. If S or R see only part of the respective
tp@0:         % edge then the edge-to-edge contribution should be "windowed" too.
tp@0:         [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:         
tp@0:     else
tp@0:         redge2_re1 = rE1_r1 + BIGEDGESTEPMATRIX(:,2)*( rE2_r1 - rE1_r1 );
tp@0:         redge1_re2 = rE1_r2 + BIGEDGESTEPMATRIX(:,1)*( rE2_r2 - rE1_r2 );                
tp@0:     end
tp@0: else
tp@0:     
tp@0:     if edgesareinplane == 0
tp@0:         error('ERROR: r and theta calculation for non-in-plane edges not implemented yet, for swapbigmatrix=1!')
tp@0: 
tp@0:     else
tp@0:         redge2_re1 = rE1_r1 + BIGE(:,2)*( rE2_r1 - rE1_r1 );
tp@0:         redge1_re2 = rE1_r2 + BIGE(:,1)*( rE2_r2 - rE1_r2 );    
tp@0:     end
tp@0: end
tp@0: if edgesareinplane == 0
tp@0: %	PS 20120415: CHECK - is it correct or not for non-in-plane edges?
tp@0: % % %     error('ERROR: r and theta calculation for non-in-plane edges not implemented yet!')
tp@0: % % %     % We need to implement something like:
tp@0: % % % %     theta = acos(r./E2Edist);
tp@0: % % %     thetaedge1_re2 = thetaE1_r2 + BIGEDGESTEPMATRIX(:,1)*( thetaE2_r2 - thetaE1_r2 );
tp@0: % % %     thetaedge2_re1 = thetaE1_r1 + BIGEDGESTEPMATRIX(:,2)*( thetaE2_r1 - thetaE1_r1 );
tp@0: else
tp@0:     thetaedge1_re2 = thetaE1_r2;
tp@0:     thetaedge2_re1 = thetaE1_r1;
tp@0: end
tp@0: 
tp@0: if method == 'n'
tp@0: 
tp@0: 	%-----------------------------------------------------------------------
tp@0: 	% In directivity function 2, DF2, we need the quantity
tp@0: 	%
tp@0: 	%   (  ( ze2_re2 - ze1_re2 ) * ( ze2_re2 - zr_re2 ) + n*l )/re1_re2/rr
tp@0: 
tp@0:     % First, E2E is relative to edge 2
tp@0:     E2Edist = sqrt( (B5 - zedge1_re2).^2 + redge1_re2.^2 );
tp@0:         
tp@0: 	% B2 will get the coshnyeta values for DF2
tp@0: 	B2 =  ( ( B5 - zedge1_re2 ).*( B5 - zR ) + E2Edist.*E2Rdist )./redge1_re2/rR;
tp@0: 	B2 = ( sqrt( B2.^2-1) + B2 ).^nyveclist(2);
tp@0: 	B2 = real( B2 + 1./B2)/2;
tp@0:     
tp@0:     clear redge1_re2 zedge1_re2
tp@0:     
tp@0: 	%-----------------------------------------------------------------------
tp@0: 	% In directivity function 1, DF1, we need the quantity
tp@0: 	%
tp@0: 	%   (  ( ze1_re1  -zs_re1 )*( ze1_re1 - ze2_re1 ) + n*m )/rs/re2_re1
tp@0: 
tp@0:     % Now, we need E2E relative to edge 1
tp@0:     E2Edist = sqrt( (B4 - zedge2_re1).^2 + redge2_re1.^2 );
tp@0: 
tp@0: 	% B1 will get the coshnyeta values for DF1
tp@0: 	
tp@0:  	B1 = (( B4 - zS ).*(  B4 - zedge2_re1 ) + S2Edist.*E2Edist)/rS./redge2_re1;
tp@0:  	B1 = ( real(sqrt( B1.^2-1)) + B1 ).^nyveclist(1);
tp@0:  	B1 = ( B1 + 1./B1)/2;
tp@0: 
tp@0:     clear zedge2_re1 redge2_re1
tp@0: 
tp@0: else
tp@0: 	B1 = 1;   B2 = 1;
tp@0: end
tp@0: 	
tp@0: % Calculate DF1.*DF2
tp@0: 
tp@0: if pathalongplane == 0
tp@0: 
tp@0:     % Note that here we use E2E re. edge 1!
tp@0:     B3 = multfac*nyveclist(1)*nyveclist(2)/16/pi^2*dzvec(1)*dzvec(2)*...
tp@0: 	              ( sin(nyveclist(1)*(pi + thetaS + thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi + thetaS + thetaedge2_re1   ))) + ...
tp@0: 				    sin(nyveclist(1)*(pi + thetaS - thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi + thetaS - thetaedge2_re1   ))) + ...
tp@0: 				    sin(nyveclist(1)*(pi - thetaS + thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi - thetaS + thetaedge2_re1   ))) + ...
tp@0: 					sin(nyveclist(1)*(pi - thetaS - thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi - thetaS - thetaedge2_re1   ))))./S2Edist./E2Edist;
tp@0: 
tp@0:     B3 = 	  B3.*( sin(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))) + ...
tp@0: 				    sin(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))) + ...
tp@0: 				    sin(nyveclist(2)*(pi - thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi - thetaedge1_re2 + thetaR))) + ...
tp@0: 					sin(nyveclist(2)*(pi - thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi - thetaedge1_re2 - thetaR))))./E2Rdist;
tp@0:                 
tp@0: else
tp@0:     
tp@0:     if thetaS == 0
tp@0:         B3 = multfac*nyveclist(1)*nyveclist(2)/2/pi^2*dzvec(1)*dzvec(2)*...
tp@0: 		              ( sin(nyveclist(1)*(pi + thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi + thetaedge2_re1   ))))./S2Edist./E2Edist;
tp@0:     else
tp@0:         
tp@0:         B3 = multfac*nyveclist(1)*nyveclist(2)/4/pi^2*dzvec(1)*dzvec(2)*...
tp@0: 		              ( sin(nyveclist(1)*(pi + thetaS + thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi + thetaS + thetaedge2_re1   ))) + ...
tp@0: 					    sin(nyveclist(1)*(pi - thetaS + thetaedge2_re1  ))./(B1- cos(nyveclist(1)*(pi - thetaS + thetaedge2_re1   ))))./S2Edist./E2Edist;
tp@0:     end    
tp@0: 
tp@0:     
tp@0:     B3 = 	      B3.*( sin(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 + thetaR))) + ...
tp@0: 					    sin(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))./(B2-cos(nyveclist(2)*(pi + thetaedge1_re2 - thetaR))))./E2Rdist;                    
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------------------------
tp@0: % Determine in which sample slots the amplitude contributions should
tp@0: % be added, based on the total distance.
tp@0: 
tp@0: % B2 is number of sample slots as non-integers
tp@0: 
tp@0: B2 = ( S2Edist + E2Edist + E2Rdist - R_irstart)/(cair/fs)+1;
tp@0: clear S2Edist E2Edist E2Rdist
tp@0: 
tp@0: % B1 is the sampleslot numbers in integer number
tp@0: B1 = floor(B2);
tp@0: 
tp@0: % B2 is a value between 0 and 1 stating how much of dh that should be added
tp@0: % to the first sample slot, i.e. the one given on B1. The rest of dh should
tp@0: % be added to the following sample slot.
tp@0: B2 = B2 - B1;
tp@0: 
tp@0: % New approach: simply state that the amplitudevalue, mult. by 1-B2, should
tp@0: % be placed in the integer slot given in B1. In addition, the
tp@0: % amplitudevalues, mult by B2, should be placed in the integer slot given
tp@0: % by B1+1.
tp@0: 
tp@0: B1 = [B1;B1+1];
tp@0: B3 = [B3.*(1-B2);B3.*B2];
tp@0: 
tp@0: ir = zeros(max(B1),1 );
tp@0: 
tp@0: [B1,sortvec] = sort(B1);
tp@0: B3 = cumsum(B3(sortvec));
tp@0: ivstep = (find(diff([(B1);(B1(end))+1])));
tp@0: ir(B1(ivstep)) = diff([0;B3(ivstep)]);
tp@0: 
tp@0: % Here we could have the possibility to save space by cutting out initial zeros and
tp@0: % set the value ninit correspondingly to the number of removed zeros. However, the rest of
tp@0: % the EDB1 functions have not implented support for this.
tp@0: 
tp@0: ninit = 0;