tp@0: function [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(ISlist,R,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec)
tp@0: % EDB1chkISvisible - Checks if paths from a set of IS to a set of R pass through their refl. planes. 
tp@0: % EDB1chkISvisible checks if the paths between a number of IS and a single R, or a number of R, 
tp@0: % pass through their respective reflecting planes.
tp@0: %
tp@0: % Input parameters:
tp@0: %   ISlist      Matrix, [nIS,3], of the nIS image source coordinates
tp@0: %   R			Receiver coordinates, [1,3] or [nIS,3]
tp@0: %   planeeqs_lastvalue  List [nIS,1] of the fourth value of the plane
tp@0: %               equations.
tp@0: %   planenvecs,minvals, maxvals, planecorners, corners, ncornersperplanevec
tp@0: %               Data that should have been taken from the corresponding
tp@0: %               variables in the eddatafile,see EDB1edgeo for more information.
tp@0: %               NB!! The matrices planeeqs, minvals, maxvals
tp@0: %               have been rearranged so that they have nIS rows, and each
tp@0: %               row contain the data for one specific plane: the reflecting
tp@0: %               plane for the IS.
tp@0: %
tp@0: % Output parameters:
tp@0: %	 hitplanes		List, [nhits,1], of the planes that have been hit. The values
tp@0: %	                in this list refer to the index numbers of the input list ISlist
tp@0: %    hitpoints      List, [nhits,3] of the hitpoint coordinates, for the planes that have been hit
tp@0: %    edgehits       List, [nedgehits,1] of the planes that were hit right at an edge.
tp@0: %                   These combinations were marked as hit in the list hitplanes, 
tp@0: %                   but the extra information in edgehits can possibly be used.
tp@0: %    edgehitpoints  List, [nedgehits,3] of the hitpoint coordinates for the
tp@0: %                   planes that were hit at an edge.
tp@0: %    cornerhits     List, [ncornerhits,1] of the planes that were hit right at a corner.
tp@0: %                   These combinations were marked as hit in the list hitplanes, 
tp@0: %                   but the extra information in edgehits can possibly be used.
tp@0: %    cornerhitpoints  List, [ncornerhits,3] of the hitpoint coordinates for the
tp@0: %                   planes that were hit at a corner.
tp@0: %
tp@0: % Uses the function EDB1poinpla.
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) 20050116
tp@0: % 
tp@0: % [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(ISlist,R,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec)
tp@0: 
tp@0: nsou = size(ISlist,1);
tp@0: nplanes = size(planeeqs_lastvalue,1);
tp@0: 
tp@0: if size(planeeqs_lastvalue,2) > 1
tp@0:     error(['Fix the call to EDB1chkISvisible!!!'])    
tp@0: end
tp@0: 
tp@0: [nrec,slask] = size(R);
tp@0: if nrec == 1
tp@0:     onesvec = uint8(ones(nplanes,1));
tp@0:     R = R(onesvec,:);
tp@0:     clear onesvec    
tp@0: end
tp@0: 
tp@0: planenvecs = planenvecs.';
tp@0: 
tp@0: %----------------------------------------------------------------
tp@0: % First we check if the IS and R are on the same side of the
tp@0: % respective planes. Then the path can not pass through the plane.
tp@0: 
tp@0: tempmatrix = corners(planecorners(:,1),:);
tp@0: sseesplanes = sum(  ((ISlist - tempmatrix).').*planenvecs   ).';
tp@0: rseesplanes = sum(  ((R - tempmatrix).').*planenvecs   ).';
tp@0: iv1 = uint32(find( (sseesplanes>= (-eps*30) )~=(rseesplanes>= (-eps*30)  ) ));
tp@0: clear sseesplanes rseesplanes
tp@0: 
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: % Next tests (if there were any planes that had IS and R at different sides)
tp@0: % 1. Are the vectors IS -> R parallel to the planes?
tp@0: % 2. Are the hitpoints inside the planes?
tp@0: 
tp@0: if ~isempty(iv1)
tp@0: 	npossible = length(iv1);
tp@0:     % Below tempmatrix is the direction vector
tp@0: 	tempmatrix = R - ISlist;
tp@0:     clear R
tp@0:     dirveclengths = sqrt(sum(tempmatrix.'.^2)).' + eps*100;
tp@0:     tempmatrix = tempmatrix./dirveclengths(:,ones(1,3));
tp@0: 
tp@0: 	% If test == 0, then the vector S -> R runs along the
tp@0: 	% plane.
tp@0: 	
tp@0: 	iv1 = iv1(find(sum(   planenvecs(:,iv1).*(tempmatrix(iv1,:).')   ).'~=0));
tp@0:     
tp@0: 	if ~isempty(iv1)
tp@0: 
tp@0: 		% The last test is if the hitpoint is inside the plane
tp@0: 	
tp@0:         % Step 1: calculate the hit point using u = the line parameter (with
tp@0:         % unit meters) from the source towards the receiver.
tp@0: 		udir = ( planeeqs_lastvalue(iv1) - sum( planenvecs(:,iv1).*(ISlist(iv1,:).') ).' );        
tp@0:         clear planeeqs_lastvalue
tp@0:         udir = udir./( sum(   planenvecs(:,iv1).*(tempmatrix(iv1,:).')   ).');
tp@0:         
tp@0:         % Step 2: check that the hit point is at a shorter distance than
tp@0:         % the receiver point, and that the hit point is not in the opposite
tp@0:         % direction than the receiver.
tp@0: 
tp@0:         iv2 = find( udir<(dirveclengths(iv1)*1.001) & udir>=0 );
tp@0:         clear dirveclengths
tp@0:         if length(iv2) < length(iv1)
tp@0:             iv1 = iv1(iv2);
tp@0:             udir = udir(iv2);
tp@0:             clear iv2
tp@0:         end
tp@0:         
tp@0:         if ~isempty(iv1)
tp@0:             % Step 3: calculate the actual xyz coordinates of the hit points
tp@0:             % Now tempmatrix gets the values of the hit points
tp@0:             tempmatrix = ISlist(iv1,:) + udir(:,ones(1,3)).*tempmatrix(iv1,:);       
tp@0:         
tp@0:             clear ISlist udir
tp@0: 
tp@0:             [hitvec,edgehitvec,cornerhitvec] = EDB1poinpla(tempmatrix,iv1,minvals,maxvals,planecorners,corners,ncornersperplanevec,planenvecs.');
tp@0:         
tp@0:     	   	hitplanes = [];
tp@0:             hitpoints = [];
tp@0:             edgehits = [];
tp@0:          	edgehitpoints = [];
tp@0:             cornerhits = [];
tp@0:             cornerhitpoints = [];
tp@0:             if any(any(hitvec)) ~=0
tp@0:                 ivhit = find(hitvec==1);
tp@0: 			    hitplanes = iv1( ivhit );  
tp@0:              hitpoints = tempmatrix(ivhit,:);
tp@0:             end
tp@0:             if ~isempty(edgehitvec)
tp@0:                 ivhit = find(edgehitvec==1);
tp@0:                 edgehits      = iv1( ivhit );
tp@0:                 edgehitpoints = tempmatrix(ivhit,:); 
tp@0:             end
tp@0:             if ~isempty(cornerhitvec)
tp@0:                 ivhit = find(cornerhitvec==1);
tp@0:                 cornerhits  = iv1( ivhit );
tp@0:                 cornerhitpoints = tempmatrix(ivhit,:);    
tp@0:             end
tp@0:         else
tp@0:     		hitplanes = [];
tp@0:             hitpoints = [];
tp@0:             edgehits = [];
tp@0:             edgehitpoints = [];
tp@0:             cornerhits = [];
tp@0:             cornerhitpoints = [];        
tp@0:         end
tp@0: 	else
tp@0: 		hitplanes = [];
tp@0:         hitpoints = [];
tp@0:         edgehits = [];
tp@0:         edgehitpoints = [];
tp@0:         cornerhits = [];
tp@0:         cornerhitpoints = [];
tp@0: 	end
tp@0: else    
tp@0: 	hitplanes = [];   
tp@0:     hitpoints = [];
tp@0:     edgehits = [];
tp@0:     edgehitpoints = [];
tp@0:     cornerhits = [];
tp@0:     cornerhitpoints = [];
tp@0: end