tp@0: function ISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfromS,vispartedgesfromS,nedgesubs)
tp@0: % EDB1findISEStree - Constructs a list ("an ISES tree") of all possible specular-diffraction combinations.
tp@0: % EDB1findISEStree constructs a list ("an ISES tree") of all possible specular-diffraction combinations that
tp@0: % are possible for the given source. This list is based on the visibility
tp@0: % matrices from source to planes/edges and between planes/edges. Visibility
tp@0: % checks are employed wherever possibble. NB! Visibility checks that depend
tp@0: % on the receiver position are not used at this stage.
tp@0: % Planes are numbered 1-nplanes, and edges have numbers
tp@0: % [nplanes+1:nplanes+nedges].
tp@0: %
tp@0: % Input parameters:
tp@0: % 	eddatafile                  See a description in EDB1edgeo, EDB1srgeo and EDB1mainISES
tp@0: % 	S
tp@0: % 	isou
tp@0: % 	specorder
tp@0: % 	difforder
tp@0: % 	visplanesfromS
tp@0: % 	vispartedgesfromS
tp@0: % 	nedgesubs
tp@0: %
tp@0: % GLOBAL parameters
tp@0: %   SHOWTEXT JJ JJnumbofchars   See EDB1mainISES
tp@0: %   SUPPRESSFILES               If this global parameter has the value 1
tp@0: %                               then the results from this function will be returned in a struct
tp@0: %                               rather than in a file.
tp@0: %
tp@0: % Output parameters:
tp@0: %   ISEStreefile                The name of the output file which contains
tp@0: %                               the output data below.
tp@0: %
tp@0: % Global output data:
tp@0: % 	POTENTIALISES               A matrix, [npossiblecombs,specorder], of all the 
tp@0: % 	                            possible specular-diffraction combos that are possible for the
tp@0: % 	                            given source. The combos are denoted by
tp@0: % 	                            plane numbers and edge numbers, but edge
tp@0: % 	                            numbers have a constant number (=nplanes)
tp@0: % 	                            added to them, i.e., if there are 20 planes
tp@0: % 	                            in the model, edge number 1 will be denoted
tp@0: % 	                            by number 21.
tp@0: % 	ORIGINSFROM                 A list, [npossiblecombs,1], where the value
tp@0: % 	                            in row N states which other row in ISCOORDS that
tp@0: % 	                            the image source in row N originates from.
tp@0: % 	ISCOORDS                    A matrix, [npossiblecombs,3], containing
tp@0: % 	                            the image source coordinates, where this is applicable
tp@0: %                               i.e., where the combo starts with a specular reflection. 
tp@0: % 	ISESVISIBILITY              A list, [npossiblecombs,1], containing the visibility of
tp@0: %                               the first edge in a sequence, seen from the source. 
tp@0: % 	IVNSPECMATRIX               A matrix, [nspeccombs,specorder], where each column contains
tp@0: %                               the row numbers in POTENTIALISES that contain combos with
tp@0: %                               1 or 2 or... or "specorder" purely specular reflections.
tp@0: % 	REFLORDER                   A list, [npossiblecombs,1], containing the
tp@0: % 	                            order of reflection for each row of POTENTIALISES
tp@0: % 	IVNDIFFMATRIX               A matrix, [ndiffcombs,specorder], where each column contains
tp@0: %                               the row numbers in POTENTIALISES that contain combos with
tp@0: %                               1 or 2 or... or "specorder" diffractions.
tp@0: % Output data in the ISEStreefile:
tp@0: % 	lengthNspecmatrix           A list, [1,specorder], with the number of
tp@0: % 	                            entries in each column of IVNSPECMATRIX.
tp@0: % 	lengthNdiffmatrix           A list, [1,specorder], with the number of
tp@0: % 	                            entries in each column of IVNDIFFMATRIX.
tp@0: % 	singlediffcol               A list, [nfirstorderdiff,1], containing the column number
tp@0: %                               that the first-order diffracted component can be found in (in
tp@0: % 	                            POTENTIALISES).
tp@0: % 	startindicessinglediff      A list, [specorder,1], where the value in position N contains
tp@0: %                               which row in IVNDIFFMATRIX(:,1) that has the first combination
tp@0: %                               of order N and with one diffraction component.
tp@0: % 	endindicessinglediff        A list, [specorder,1], where the value in position N contains
tp@0: %                               which row in IVNDIFFMATRIX(:,1) that has the last combination
tp@0: %                               of order N and with one diffraction 
tp@0: % 	ndecimaldivider             See below.
tp@0: % 	PointertoIRcombs            A sparse list which for position N contains a row number
tp@0: %                               (in POTENTIALISES) that has a combination
tp@0: %                               that ends with a specular reflection in
tp@0: %                               plane N, after a diffraction. A row with one of the specular
tp@0: %                               reflections M,N after a diffraction can be
tp@0: %                               found in position M*ndecimaldivider+N etc
tp@0: %                               for higher orders.
tp@0: % 	IRoriginsfrom               A list, [npossiblecombs,1], where the value in position N
tp@0: %                               states which other row number (in POTENTIALISES) that the
tp@0: %                               image receiver in row N origins from.
tp@0: %
tp@0: % Uses functions EDB1strpend EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths
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) 20100210
tp@0: %
tp@0: % ISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms,vispartedgesfroms,nedgesubs);
tp@0: 
tp@0: global SHOWTEXT JJ JJnumbofchars
tp@0: global POTENTIALISES ISCOORDS IVNDIFFMATRIX
tp@0: global IVNSPECMATRIX ORIGINSFROM ISESVISIBILITY REFLORDER SUPPRESSFILES
tp@0: 
tp@0: GEOMACC = 1e-10;
tp@0: 
tp@0: eval(['load ',eddatafile])
tp@0: clear cornerinfrontofplane
tp@0: 
tp@0: [filepath,filestem,fileext] = fileparts(eddatafile);
tp@0: ISEStreefile = [[filepath,filesep],EDB1strpend(filestem,'_eddata'),'_',int2str(isou),'_ISEStree.mat'];
tp@0: 
tp@0: ncorners = size(corners,1);
tp@0: nplanes = length(planeisthin);
tp@0: nedges = length(closwedangvec);
tp@0: nhighval = nplanes + nedges;
tp@0: 
tp@0: onesvec1 = ones(1,nedgesubs);
tp@0: obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0: onesvec3 = ones(1,3);
tp@0: 
tp@0: doconetrace = 0;
tp@0: if SHOWTEXT >= 3
tp@0:     disp('WARNING: doconetrace is switched off')    
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: % We modify edgeseesplane because in EDB1edgeo, edgeseesplane was set to 1
tp@0: % for totabs planes (this was necessary for EDB1srgeo).
tp@0: 
tp@0: % Set all edges belonging to planes with reflfactors = 0 to edgeseesplane = 0.
tp@0: 
tp@0: listofabsorbingplanes = find(reflfactors == 0);
tp@0: 
tp@0: if ~isempty(listofabsorbingplanes)
tp@0: 	edgeseesplane(listofabsorbingplanes,:) = zeros(length(listofabsorbingplanes),nedges);
tp@0: end
tp@0: 
tp@0: % If we should do cone tracing, then we need plane midpoints and plane
tp@0: % sizes.
tp@0: 
tp@0: if doconetrace == 1
tp@0:     planemidpoints = zeros(nplanes,3);
tp@0:     maxdisttocorner = zeros(nplanes,1);
tp@0: 
tp@0:     iv4planes = find(ncornersperplanevec==4);
tp@0:     if any(iv4planes)
tp@0:         xvalues = corners(planecorners(iv4planes,1:4).',1);
tp@0:         xvalues = reshape(xvalues,4,length(xvalues)/4);
tp@0:         yvalues = corners(planecorners(iv4planes,1:4).',2);
tp@0:         yvalues = reshape(yvalues,4,length(yvalues)/4);
tp@0:         zvalues = corners(planecorners(iv4planes,1:4).',3);
tp@0:         zvalues = reshape(zvalues,4,length(zvalues)/4);
tp@0:         planemidpoints(iv4planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0:         
tp@0:         dist1 = sqrt(sum((corners(planecorners(iv4planes,1),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0:         dist2 = sqrt(sum((corners(planecorners(iv4planes,2),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0:         dist3 = sqrt(sum((corners(planecorners(iv4planes,3),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0:         dist4 = sqrt(sum((corners(planecorners(iv4planes,4),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0:         maxdisttocorner(iv4planes) = max([dist1 dist2 dist3 dist4].');
tp@0:         
tp@0:     end
tp@0: 
tp@0:     iv3planes = find(ncornersperplanevec==3);
tp@0:     if any(iv3planes)
tp@0:         xvalues = corners(planecorners(iv3planes,1:3).',1);
tp@0:         xvalues = reshape(xvalues,3,length(xvalues)/3);
tp@0:         yvalues = corners(planecorners(iv3planes,1:3).',2);
tp@0:         yvalues = reshape(yvalues,3,length(yvalues)/3);
tp@0:         zvalues = corners(planecorners(iv3planes,1:3).',3);
tp@0:         zvalues = reshape(zvalues,3,length(zvalues)/3);
tp@0:         planemidpoints(iv3planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0: 
tp@0:         dist1 = sqrt(sum((corners(planecorners(iv3planes,1),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0:         dist2 = sqrt(sum((corners(planecorners(iv3planes,2),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0:         dist3 = sqrt(sum((corners(planecorners(iv3planes,3),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0:         maxdisttocorner(iv3planes) = max([dist1 dist2 dist3].');
tp@0:     
tp@0:     end
tp@0: 
tp@0:     iv5planes = find(ncornersperplanevec==5);
tp@0:     if any(iv5planes)
tp@0:         xvalues = corners(planecorners(iv5planes,1:5).',1);
tp@0:         xvalues = reshape(xvalues,5,length(xvalues)/5);
tp@0:         yvalues = corners(planecorners(iv5planes,1:5).',2);
tp@0:         yvalues = reshape(yvalues,5,length(yvalues)/5);
tp@0:         zvalues = corners(planecorners(iv5planes,1:5).',3);
tp@0:         zvalues = reshape(zvalues,5,length(zvalues)/5);
tp@0:         planemidpoints(iv5planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0:     
tp@0:         dist1 = sqrt(sum((corners(planecorners(iv5planes,1),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0:         dist2 = sqrt(sum((corners(planecorners(iv5planes,2),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0:         dist3 = sqrt(sum((corners(planecorners(iv5planes,3),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0:         dist4 = sqrt(sum((corners(planecorners(iv5planes,4),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0:         dist5 = sqrt(sum((corners(planecorners(iv5planes,5),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0:         maxdisttocorner(iv5planes) = max([dist1 dist2 dist3 dist4 dist5].');
tp@0: 
tp@0:     end
tp@0: 
tp@0:     iv6planes = find(ncornersperplanevec==6);
tp@0:     if any(iv6planes)
tp@0:         xvalues = corners(planecorners(iv6planes,1:6).',1);
tp@0:         xvalues = reshape(xvalues,6,length(xvalues)/6);
tp@0:         yvalues = corners(planecorners(iv6planes,1:6).',2);
tp@0:         yvalues = reshape(yvalues,6,length(yvalues)/6);
tp@0:         zvalues = corners(planecorners(iv6planes,1:6).',3);
tp@0:         zvalues = reshape(zvalues,6,length(zvalues)/6);
tp@0:         planemidpoints(iv6planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0: 
tp@0:         dist1 = sqrt(sum((corners(planecorners(iv6planes,1),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         dist2 = sqrt(sum((corners(planecorners(iv6planes,2),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         dist3 = sqrt(sum((corners(planecorners(iv6planes,3),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         dist4 = sqrt(sum((corners(planecorners(iv6planes,4),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         dist5 = sqrt(sum((corners(planecorners(iv6planes,5),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         dist6 = sqrt(sum((corners(planecorners(iv6planes,6),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0:         maxdisttocorner(iv6planes) = max([dist1 dist2 dist3 dist4 dist5 dist6].');
tp@0:     
tp@0:     end
tp@0: 
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: % Combine planeseesplane with edgeseesplane
tp@0: % and visplanesfroms with vispartedgesfroms
tp@0: % so that edges appear as planes
tp@0: 
tp@0: if difforder > 0
tp@0:     visPLANESfroms = [visplanesfromS;2*sign(double(vispartedgesfromS))];
tp@0:     % source sees planes that have visplanesfroms:
tp@0:     %   2 => in front of reflective planes
tp@0:     %   4 => inside plane which is reflective
tp@0:     %   old version:  sousees = (visPLANESfroms==1 | visPLANESfroms==-2);
tp@0:     sousees = (visPLANESfroms==2 | visPLANESfroms==4);
tp@0: 
tp@0:     if exist('edgeseespartialedge') ~= 1
tp@0:         edgeseesedge = int8(zeros(nedges,nedges));    
tp@0:     else
tp@0:         if ~isempty(edgeseespartialedge)
tp@0: 
tp@0:             edgeseesedge = int8(full(abs(double(edgeseespartialedge))>0));
tp@0:             
tp@0:         else
tp@0:             edgeseesedge = int8(zeros(nedges,nedges));    
tp@0:         end
tp@0:     end
tp@0: 
tp@0:     PLANEseesPLANE = [[planeseesplane edgeseesplane];[edgeseesplane.' edgeseesedge]];
tp@0:     
tp@0:     clear edgeseesplane edgeseesedge
tp@0: else
tp@0:     visPLANESfroms = visplanesfromS;
tp@0:     PLANEseesPLANE = planeseesplane;
tp@0:     sousees = (visPLANESfroms==2 | visPLANESfroms==4);
tp@0: 
tp@0:     clear visplanesfromS
tp@0: end
tp@0: 
tp@0: startindices = zeros(specorder,1);
tp@0: endindices = zeros(specorder,1);
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: % Construct a list of which planes a plane sees so that a search isn't
tp@0: % needed later.
tp@0: 
tp@0: nPLANES = size(PLANEseesPLANE,1);
tp@0: listofvisPLANES = zeros(nPLANES,nPLANES-1);    
tp@0: nvisPLANES = zeros(nPLANES,1);
tp@0: for ii = 1:nPLANES
tp@0:     visPLANES = find(PLANEseesPLANE(ii,:)==1);
tp@0:     nvisPLANES(ii) = length(visPLANES);
tp@0:     listofvisPLANES(ii,1:nvisPLANES(ii)) = visPLANES;
tp@0: end
tp@0: 
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %
tp@0: %         First order
tp@0: %
tp@0: %------------------------------------------------------------------
tp@0: 
tp@0: startindices(1) = 1;
tp@0: 
tp@0: possiblefirsts  =  find( sousees );
tp@0: npossiblefirsts = length(possiblefirsts);
tp@0: 
tp@0: if nhighval < 255
tp@0:     POTENTIALISES = uint8( [[possiblefirsts ] zeros(npossiblefirsts,specorder-1)] );
tp@0: else
tp@0:     POTENTIALISES = uint16( [[possiblefirsts ] zeros(npossiblefirsts,specorder-1)] );
tp@0: end
tp@0: 
tp@0: ORIGINSFROM = uint32(zeros(npossiblefirsts,1));
tp@0: if nedgesubs < 8 
tp@0:     ISESVISIBILITY = uint8(zeros(npossiblefirsts,1));
tp@0: elseif nedgesubs < 16
tp@0:     ISESVISIBILITY = uint16(zeros(npossiblefirsts,1));
tp@0: else    
tp@0:     ISESVISIBILITY = uint32(zeros(npossiblefirsts,1));
tp@0: end
tp@0: iv = find(possiblefirsts>nplanes);
tp@0: if ~isempty(iv)
tp@0:     ISESVISIBILITY(iv) = vispartedgesfromS(possiblefirsts(iv)-nplanes);    
tp@0: end
tp@0: 
tp@0: endindices(1) = npossiblefirsts;
tp@0: 
tp@0: % Compute the IS coords
tp@0: 
tp@0: ivdiff = [startindices(1):endindices(1)];
tp@0: ivspec = find(POTENTIALISES(ivdiff,1)<=nplanes);
tp@0: ivdiff(ivspec) = [];
tp@0: 
tp@0: ISCOORDS = zeros(npossiblefirsts,3);
tp@0: ISCOORDS(ivspec,:) = EDB1findis(S,POTENTIALISES(ivspec,1),planeeqs,1,onesvec3);
tp@0: 
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %
tp@0: % Second order
tp@0: %
tp@0: %---------------------------------------------------------------------------------------------
tp@0: 
tp@0: startindices(2) = startindices(1) + length(possiblefirsts);
tp@0: 
tp@0: if specorder > 1
tp@0:     if SHOWTEXT >= 3
tp@0:         disp(['     Order number 2'])
tp@0:     end
tp@0:     
tp@0:     if nhighval < 255
tp@0:         startplanes = uint8(possiblefirsts);
tp@0:     else
tp@0:         startplanes = uint16(possiblefirsts);
tp@0:     end
tp@0:         
tp@0: 	addtocomb = startplanes;
tp@0: 	nadds = size(addtocomb,1);
tp@0:     if nadds > 0
tp@0: 
tp@0:        	maxnumberofvispla = max(sum(PLANEseesPLANE(:,startplanes)==1));
tp@0: 	
tp@0:         addtocomb = (reshape(addtocomb(:,ones(1,maxnumberofvispla)).',length(addtocomb)*maxnumberofvispla,1));
tp@0: 		addtocomb = [addtocomb zeros(size(addtocomb))];
tp@0:         
tp@0:         addtoISESVISIBILITY = reshape(ISESVISIBILITY(:,ones(1,maxnumberofvispla)).',nadds*maxnumberofvispla,1);
tp@0:         
tp@0:         startpos = [[1:maxnumberofvispla:length(addtocomb)]  length(addtocomb)+1  ].';
tp@0: 		
tp@0: 		for ii = 1:length(startpos)-1
tp@0:             if nvisPLANES(addtocomb(startpos(ii))) > 0
tp@0:     			possibleplanes = listofvisPLANES(addtocomb(startpos(ii)),1:nvisPLANES(addtocomb(startpos(ii)))).';
tp@0: 			    nposs = length(possibleplanes);
tp@0: 			    addtocomb( startpos(ii):startpos(ii)+nposs-1,2) = possibleplanes;
tp@0:             end
tp@0:             
tp@0: 		end
tp@0: 	
tp@0:         addtooriginsfrom = uint32([1:startindices(2)-1].');
tp@0: 		addtooriginsfrom = addtooriginsfrom(:,ones(1,maxnumberofvispla));
tp@0: 		addtooriginsfrom = reshape(addtooriginsfrom.',maxnumberofvispla*(startindices(2)-1),1);
tp@0: 		
tp@0: 		cutout = find(addtocomb(:,2)==0);
tp@0: 		addtocomb(cutout,:) = [];
tp@0: 		addtooriginsfrom(cutout) = [];
tp@0:         addtoISESVISIBILITY(cutout) = [];
tp@0:         clear cutout
tp@0:         nadds = size(addtocomb,1);
tp@0: 		    
tp@0:         if nadds > 0
tp@0: 
tp@0:             %--------------------------------------------------------------------------
tp@0: 			% Try to get rid of some combinations that we know can not be propagated
tp@0:                         
tp@0: 			% Find those combinations of specular-specular where the IS is behind the second reflecting plane
tp@0: 			% because they can then never ever be propagated, or valid
tp@0: 			
tp@0: 			if difforder > 0
tp@0:                 ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0: 			else
tp@0:                 ivss = uint32([1:nadds].');
tp@0: 			end
tp@0: 
tp@0: 			imsoucheck = dot((ISCOORDS(addtooriginsfrom(ivss),1:3) - corners(planecorners(addtocomb(ivss,2),1),:)).',planenvecs(addtocomb(ivss,2),:).').';
tp@0: 			imsounotOK = find(imsoucheck < GEOMACC);
tp@0: 			addtocomb(ivss(imsounotOK),:) = [];
tp@0: 			addtooriginsfrom(ivss(imsounotOK),:) = [];
tp@0:             addtoISESVISIBILITY(ivss(imsounotOK)) = [];
tp@0: 
tp@0:             nadds = size(addtocomb,1);
tp@0: 
tp@0:             if nadds > 0 & doconetrace == 1
tp@0: 
tp@0:                 if difforder > 0
tp@0:                     ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0: 				else
tp@0:                     ivss = uint32([1:nadds].');
tp@0: 				end
tp@0:                 
tp@0:                 beamdirection1 = planemidpoints(addtocomb(ivss,1),:) - ISCOORDS(addtocomb(ivss,1),:);
tp@0:                 beamlength = sqrt( sum( beamdirection1.'.^2 ) ).';
tp@0:                 beamdirection1  = beamdirection1./beamlength(:,ones(1,3));                
tp@0:                 maxprojradius = maxdisttocorner(addtocomb(ivss,1));
tp@0:                 % The maxdisttocorner below should be multiplied
tp@0:                 % by sine of the angle between the plane normal vector and the
tp@0:                 % beam direction.
tp@0:                 iv = find(beamlength<maxprojradius);
tp@0:                 if ~isempty(iv)
tp@0:                     beamlength(iv) = beamlength(iv)*0+ eps*1e6;    
tp@0:                 end
tp@0:                 beamangle1 = atan(maxprojradius./beamlength);
tp@0:                 
tp@0:                 % Now we calculate the beam-widths to the second reflection
tp@0:                 % plane, seen from the first-order image source. If this
tp@0:                 % beamwidth is outside the previously calculated beam of
tp@0:                 % reflection plane 1, then plane 2 can never be seen
tp@0:                 % through plane 1.
tp@0: 
tp@0:                 beamdirection2 = planemidpoints(addtocomb(ivss,2),:) - ISCOORDS(addtocomb(ivss,1),:);
tp@0:                 beamlength = sqrt( sum( beamdirection2.'.^2 ) ).';
tp@0:                 beamdirection2  = beamdirection2./beamlength(:,ones(1,3));                
tp@0:                 maxprojradius = maxdisttocorner(addtocomb(ivss,2));
tp@0:                 % The maxdisttocorner below should be multiplied
tp@0:                 % by sine of the angle between the plane normal vector and the
tp@0:                 % beam direction.
tp@0:                 iv = find(beamlength<maxprojradius);
tp@0:                 if ~isempty(iv)
tp@0:                     beamlength(iv) = beamlength(iv)*0+ eps*10;    
tp@0:                 end
tp@0:                 beamangle2 = atan(maxprojradius./beamlength);
tp@0:                 anglebetweenbeams = acos(sum( (beamdirection1.*beamdirection2).' ).');
tp@0:                 ivcutout = find( (beamangle1 + beamangle2 < anglebetweenbeams) & (beamangle1 <1.56) & (beamangle2 < 1.56));
tp@0:                 if ~isempty(ivcutout)
tp@0: 					addtocomb(ivss(ivcutout),:) = [];
tp@0: 					addtooriginsfrom(ivss(ivcutout),:) = [];
tp@0:                     addtoISESVISIBILITY(ivss(ivcutout)) = [];
tp@0:                     nadds = size(addtocomb,1);
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:             
tp@0: 			% Combinations of diffraction-specular don't need to be checked because
tp@0: 			% 'edgeseesplane' should have taken care of this.
tp@0: 			%%%%ivds = find(addtocomb(:,1)>nplanes  & addtocomb(:,2)<=nplanes);
tp@0: 			
tp@0: 			%--------------------------------------------------------------------------
tp@0: 			% Combinations of specular-diffraction: check if the IS is behind both of
tp@0: 			% the two planes that the reflecting edge connects to.
tp@0:             %
tp@0:             % Bug found 20080711 PS: The visibility test must be split into two, depending on the wedge angle!! 
tp@0:             % If the open wedge angle > pi, then it is correct to check if the IS behind both of the planes. 
tp@0:             % If the open wedge angle < pi, then it should be checked if the IS behind one of the planes!! 
tp@0:             
tp@0:             if difforder > 0 & nadds > 0
tp@0: 
tp@0:                 % First we check the sd combinations where the open wedge angle > pi
tp@0:                 % For those cases, we can exclude combinations where the IS
tp@0:                 % is behind both planes
tp@0:                 
tp@0:   				ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0:                 ivsd1 = ivsd( find(   closwedangvec( double(addtocomb(ivsd,2))-nplanes ) < pi ));
tp@0: 
tp@0:                 if ~isempty(ivsd1)
tp@0:                     imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,1)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,1)    ,:).').';
tp@0:                     imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,2)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,2)    ,:).').';
tp@0:                     GEOMACC = 1e-10;
tp@0:                     imsounotOK = find(imsoucheck1 < GEOMACC & imsoucheck2 < GEOMACC);
tp@0:                     clear imsoucheck1 imsoucheck2
tp@0: 
tp@0:                     addtocomb(ivsd1(imsounotOK),:) = [];
tp@0:                     addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0:                     addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0:                   end
tp@0:                 
tp@0:                 % Second we check the sd combinations where the open wedge angle < pi
tp@0:                 % For those cases, we can exclude combinations where the IS
tp@0:                 % is behind one of the two planes
tp@0:                 
tp@0:   				ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0:                 if ~isempty(ivsd)
tp@0:                     ivsd1 = ivsd( find(   closwedangvec( double(addtocomb(ivsd,2))-nplanes ) > pi ));
tp@0: 
tp@0:                     if ~isempty(ivsd1)
tp@0:                         imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,1)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,1)    ,:).').';
tp@0:                         imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,2)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,2))-nplanes,2)    ,:).').';
tp@0:                         GEOMACC = 1e-10;
tp@0:                         imsounotOK = find(imsoucheck1 < GEOMACC | imsoucheck2 < GEOMACC);
tp@0:                         clear imsoucheck1 imsoucheck2
tp@0: 
tp@0:                         addtocomb(ivsd1(imsounotOK),:) = [];
tp@0:                         addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0:                         addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0:                     end
tp@0:                 end
tp@0:                 
tp@0: 				nadds = size(addtocomb,1);
tp@0:                 
tp@0:                 ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0:                 
tp@0:                 if ~isempty(ivsd)
tp@0: 		
tp@0:                     %----------------------------------------------------------------------
tp@0:                     % Combinations of specular-diffraction that are not visible at
tp@0:                     % all should be removed and then they are not propagated either.
tp@0: 				
tp@0:                     if nedges < 255
tp@0:                         possibleedges = uint8(double(addtocomb(ivsd,2)) - nplanes);
tp@0:                     else
tp@0:                         possibleedges = uint16(double(addtocomb(ivsd,2)) - nplanes);
tp@0:                     end
tp@0:                      
tp@0:                     possiblecombs = addtocomb(ivsd,1);
tp@0:                     reftoISCOORDS = addtooriginsfrom(ivsd);
tp@0: 				
tp@0:                     % Expand to take the edge subdivisions into account
tp@0:                     
tp@0:                     nposs = length(ivsd);
tp@0:                     nposs = nposs*nedgesubs;        % We treat all the little edge subdivisions as separate edges
tp@0: 				
tp@0:                     expandedposscombs = possiblecombs(:,onesvec1);
tp@0:                     clear possiblecombs
tp@0:                     expandedposscombs = reshape(expandedposscombs.',nposs,1);			
tp@0:                     
tp@0:                     expandedreftoISCOORDS = reftoISCOORDS(:,onesvec1);
tp@0:                     expandedreftoISCOORDS = reshape(expandedreftoISCOORDS.',nposs,1);
tp@0:                     expandedpossedges = possibleedges(:,onesvec1);
tp@0:                     expandedpossedges = reshape(expandedpossedges.',nposs,1);
tp@0:                     expandedivsd = ivsd(:,onesvec1);
tp@0:                     expandedivsd = reshape(expandedivsd.',nposs,1);
tp@0:                     
tp@0:                     if SHOWTEXT >= 3
tp@0: 						disp(['         ',int2str(nposs),' IS+edge segments found initially '])
tp@0: 					end
tp@0:                     
tp@0:                     if nposs > 0
tp@0: 				
tp@0:                         fromcoords = ISCOORDS(expandedreftoISCOORDS,:);
tp@0:                         [tocoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0:                         clear possibleedges
tp@0:                         
tp@0:                          [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,1),4),planenvecs(expandedposscombs(:,1),:),minvals(expandedposscombs(:,1),:),...
tp@0:  						    maxvals(expandedposscombs(:,1),:),planecorners(expandedposscombs(:,1),:),corners,ncornersperplanevec(expandedposscombs(:,1)));
tp@0:                         if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0:                             disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0:                             disp('         handled correctly yet.')
tp@0:                         end
tp@0: 
tp@0:                         eval(['reflpoints',JJ(1,1),' = reflpoints;'])
tp@0: 				
tp@0:                         expandedivsd          = expandedivsd(hitplanes);
tp@0:                         expandedposscombs     = expandedposscombs(hitplanes,:);
tp@0:                         expandedreftoISCOORDS = expandedreftoISCOORDS(hitplanes);
tp@0:                         expandedpossedges     = expandedpossedges(hitplanes);
tp@0:                         edgeweightlist        = edgeweightlist(hitplanes);
tp@0:                         toedgecoords          = tocoords(hitplanes,:);
tp@0: 				
tp@0:                         nposs = length(expandedivsd);
tp@0: 				
tp@0:                     end
tp@0:                     
tp@0:                     if SHOWTEXT >= 3
tp@0: 						disp(['         ',int2str(nposs),' IS+edge segments visible '])
tp@0: 					end
tp@0: 			
tp@0:                     if obstructtestneeded & nposs > 0
tp@0:                         
tp@0:                         for jj = 1:2
tp@0:                             if nposs > 0
tp@0:                                 if jj==1
tp@0:                                     fromcoords = S;    
tp@0:                                     startplanes = [];    
tp@0:                                 else
tp@0:                                     startplanes = expandedposscombs(:,jj-1);
tp@0:                                     eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0:                                 end
tp@0:                                 if jj == 2
tp@0:                                     tocoords = toedgecoords;
tp@0:                                     endplanes = [planesatedge(expandedpossedges,1) planesatedge(expandedpossedges,2)];
tp@0:                                 else
tp@0:                                     eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])    
tp@0:                                     endplanes = expandedposscombs(:,jj);    
tp@0:                                 end
tp@0: 
tp@0:                                 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0:                                     planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0:                                 
tp@0:                                 if nobstructions > 0
tp@0:                                     expandedivsd          = expandedivsd(nonobstructedpaths);
tp@0:                                     expandedposscombs     = expandedposscombs(nonobstructedpaths,:);
tp@0:                                     expandedreftoISCOORDS = expandedreftoISCOORDS(nonobstructedpaths);
tp@0:                                     expandedpossedges     = expandedpossedges(nonobstructedpaths);
tp@0:                                     edgeweightlist        = edgeweightlist(nonobstructedpaths);
tp@0:                                     toedgecoords          = tocoords(nonobstructedpaths,:);
tp@0:                                     nposs = length(expandedivsd);
tp@0:                                     
tp@0:                                     for kk = 1:1, %norder
tp@0:                                         eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])    
tp@0:                                     end
tp@0:                                     
tp@0:                                 end
tp@0:                             end
tp@0:                  
tp@0:                         end
tp@0:                     end        
tp@0:                     
tp@0:                     if SHOWTEXT >= 3
tp@0: 					    disp(['         ',int2str(nposs),' IS+edge segments survived the obstruction test'])
tp@0:                     end
tp@0:                     
tp@0:                     % There are repetitions of each combination since each edge segment has
tp@0:                     % been treated separately. Pack them together now
tp@0: 			
tp@0:                     test = [expandedposscombs expandedpossedges];
tp@0:                     ncombs = length(expandedpossedges);
tp@0:                     dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
tp@0:                     ivremove = find(dtest==1);
tp@0:                     
tp@0:                     while ~isempty(ivremove) & ~isempty(edgeweightlist),
tp@0:                         edgeweightlist(ivremove+1) = double(edgeweightlist(ivremove+1)) + double(edgeweightlist(ivremove));
tp@0:                         edgeweightlist(ivremove) = [];
tp@0:                         expandedpossedges(ivremove) = [];
tp@0:                         expandedposscombs(ivremove,:) = [];
tp@0:                         expandedivsd(ivremove) = [];
tp@0:                         nposs = length(expandedivsd);   
tp@0:                         
tp@0:                         test = [expandedposscombs expandedpossedges];
tp@0:                         ncombs = length(expandedpossedges);
tp@0:                         dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
tp@0:                         ivremove = find(dtest==1);
tp@0: 			
tp@0:                     end
tp@0: 				       
tp@0:                     if SHOWTEXT >= 3
tp@0: 						disp(['         ',int2str(nposs),' IS+edge combinations after edge segments are packed together'])
tp@0:                     end
tp@0:                     combstoremove = setdiff(ivsd,expandedivsd);
tp@0:                     addtocomb(combstoremove,:) = [];
tp@0:                     addtooriginsfrom(combstoremove) = [];  
tp@0:                     addtoISESVISIBILITY(combstoremove) = [];
tp@0:                     nadds = length(addtooriginsfrom);
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:             
tp@0: 			% Combinations of diffraction-diffraction don't need to be checked because
tp@0: 			% 'edgeseesedge' should have taken care of this.
tp@0: 			
tp@0:             %----------------------------------------------------------------------
tp@0:             % Add the new combinations to the master list
tp@0: 		
tp@0:             if nadds > 0
tp@0:                 POTENTIALISES = [POTENTIALISES;[addtocomb zeros(nadds,specorder-2)]];
tp@0: 				ORIGINSFROM = [ORIGINSFROM;addtooriginsfrom];
tp@0: 			                
tp@0:                 if difforder > 0
tp@0:                     ivtemp = find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes);
tp@0:                     if ~isempty(ivtemp),        
tp@0:                         addtoISESVISIBILITY(ivtemp) = edgeweightlist;
tp@0:                     end
tp@0:                     ivtemp = find(addtocomb(:,1)>nplanes);
tp@0:                 end
tp@0:                 ISESVISIBILITY = [ISESVISIBILITY;addtoISESVISIBILITY];                
tp@0: 
tp@0: 				endindices(2) = length(ORIGINSFROM);
tp@0: 			
tp@0: 				% Compute the IS coords of the combinations with only specular reflections
tp@0: 				
tp@0: 				ISCOORDStoadd = zeros(nadds,3);
tp@0: 				if difforder > 0
tp@0:                     ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0: 				else
tp@0:                     ivss = uint32([1:nadds].');    
tp@0: 				end
tp@0: 				soutoreflect = ISCOORDS(ORIGINSFROM(double(ivss)+endindices(1)),1:3);
tp@0: 				ISCOORDStoadd(ivss,:) = EDB1findis(soutoreflect,POTENTIALISES(double(ivss)+endindices(1),2),planeeqs,size(soutoreflect,1),onesvec3);
tp@0:                 clear soutoreflect
tp@0: 				ISCOORDS = [ISCOORDS;ISCOORDStoadd];
tp@0: 		
tp@0:             end
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %##################################################################
tp@0: %
tp@0: % Third order and higher
tp@0: %
tp@0: %---------------------------------------------------------------------------------------------
tp@0: 
tp@0: for ordernum = 3:specorder
tp@0:     if SHOWTEXT >= 3
tp@0:         disp(['     Order number ',int2str(ordernum)])
tp@0:     end
tp@0:         
tp@0:     startindices(ordernum) = startindices(ordernum-1) + length(addtocomb);
tp@0: 
tp@0:     % The lines below could use
tp@0:     % planestoprop = unique(addtocomb(:,ordernum-1));
tp@0:     
tp@0:     planelist = sort(addtocomb(:,ordernum-1));
tp@0:     
tp@0:     if ~isempty(planelist),            
tp@0:         planestoprop = planelist([1;find(diff(double(planelist))~=0)+1]);
tp@0:     else
tp@0:         planestoprop = [];    
tp@0:     end
tp@0:     clear planelist
tp@0: 
tp@0: 	maxnumberofvispla = max(sum(PLANEseesPLANE(:,planestoprop)==1));
tp@0:     oldaddtocomb = addtocomb;
tp@0: 
tp@0: 	nadds = size(addtocomb,1);
tp@0:     if nhighval < 255
tp@0:         addtocomb = uint8(zeros(nadds*maxnumberofvispla,ordernum));
tp@0:     else
tp@0:         addtocomb = uint16(zeros(nadds*maxnumberofvispla,ordernum));
tp@0:     end
tp@0: 
tp@0:     onesvec2 = ones(1,maxnumberofvispla);
tp@0: 	for ii = 1:ordernum-1
tp@0: 		temp = oldaddtocomb(:,ii);
tp@0: 		addtocomb(:,ii) = reshape(temp(:,onesvec2).',length(temp)*maxnumberofvispla,1);
tp@0: 	end
tp@0:     nrows = size(addtocomb,1);
tp@0: 	clear oldaddtocomb temp
tp@0:     
tp@0:     if nrows > 0
tp@0:        	startpos = [[1:maxnumberofvispla:nrows]  nrows+1  ].';
tp@0:     else
tp@0:         startpos = 1;    
tp@0:     end
tp@0: 
tp@0: 	for ii = 1:length(startpos)-1
tp@0:             if nvisPLANES(addtocomb(startpos(ii),ordernum-1)) > 0
tp@0:     			possibleplanes = listofvisPLANES(addtocomb(startpos(ii),ordernum-1),1:nvisPLANES(addtocomb(startpos(ii),ordernum-1))).';
tp@0: 			    nposs = length(possibleplanes);
tp@0: 			    addtocomb( startpos(ii):startpos(ii)+nposs-1,ordernum) = possibleplanes;
tp@0:             end
tp@0: 	end
tp@0: 
tp@0:     addtooriginsfrom = uint32([1:nadds].' + startindices(ordernum-1)-1);
tp@0:     addtooriginsfrom = addtooriginsfrom(:,ones(1,maxnumberofvispla));
tp@0:     addtooriginsfrom = reshape(addtooriginsfrom.',maxnumberofvispla*nadds,1);
tp@0: 
tp@0:     addtoISESVISIBILITY = reshape(addtoISESVISIBILITY(:,ones(1,maxnumberofvispla)).',nadds*maxnumberofvispla,1);
tp@0:     
tp@0: 	clear startpos
tp@0: 	cutout = find(addtocomb(:,ordernum)==0);
tp@0: 	addtocomb(cutout,:) = [];
tp@0:     addtooriginsfrom(cutout) = [];
tp@0:     addtoISESVISIBILITY(cutout) = [];
tp@0:     nadds = size(addtocomb,1);
tp@0:     
tp@0: 	if ordernum == specorder
tp@0: 		clear cutout
tp@0: 	end
tp@0: 
tp@0:     %--------------------------------------------------------------------------
tp@0:     % Try to get rid of some combinations that we know can not be propagated
tp@0: 
tp@0:     % Find those combinations of specular-specular where the IS is behind the second reflecting plane
tp@0:     % because they can then never ever be propagated.
tp@0: 
tp@0:     if difforder > 0
tp@0:         ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0:     else
tp@0:         ivss = uint32([1:nadds].');
tp@0:     end
tp@0:     
tp@0:     imsoucheck = dot((ISCOORDS(addtooriginsfrom(ivss),1:3) - corners(planecorners(addtocomb(ivss,ordernum),1),:)).',planenvecs(addtocomb(ivss,ordernum),:).').';
tp@0: 	GEOMACC = 1e-10;
tp@0: 	imsounotOK = uint32(find(imsoucheck < GEOMACC));
tp@0:     clear imsoucheck
tp@0: 	addtocomb(ivss(imsounotOK),:) = [];
tp@0: 	addtooriginsfrom(ivss(imsounotOK),:) = [];
tp@0:     addtoISESVISIBILITY(ivss(imsounotOK)) = [];
tp@0:     clear imsounotOK
tp@0: 
tp@0:     % Addition 20100210 PS: we should check all combinations with dss in the
tp@0:     % three last orders of addtocomb: if the edge in the first 'd' belongs to the
tp@0:     % plane of the last 's' *and* the two 's' planes form a 90 degree
tp@0:     % corner, then these combos should be removed.
tp@0:     % 
tp@0:     % As a first step, we check if there are any interior-90-degree
tp@0:     % corners/edges of the model, because only then can these combos occur.
tp@0:     
tp@0:     if ordernum == 3
tp@0:        deg90edges = find( abs(closwedangvec-3*pi/2) < GEOMACC );
tp@0:        if any(deg90edges)
tp@0:            deg90planepairs = planesatedge(deg90edges,:);    
tp@0:        end
tp@0:     end    
tp@0:     
tp@0:     if difforder > 0 & any(deg90edges)
tp@0:         ivdss = uint32(find(    double(addtocomb(:,ordernum-2)>nplanes).*double(addtocomb(:,ordernum-1)<=nplanes).*double(addtocomb(:,ordernum)<=nplanes)      ));
tp@0:     else
tp@0:         ivdss = [];
tp@0:     end
tp@0: 
tp@0:     if ~isempty(ivdss)
tp@0:         
tp@0:         % First we check if the 'd' (the edge) belongs to the last 's'
tp@0:         A = addtocomb(ivdss,ordernum-2:ordernum);
tp@0: 
tp@0:         ivsubset = find( (planesatedge(A(:,1)-nplanes,1) == A(:,3)) | (planesatedge(A(:,1)-nplanes,2) == A(:,3)) );
tp@0:         ivdss = ivdss(ivsubset);
tp@0: 
tp@0:         if ~isempty(ivdss)
tp@0:             % Second, we check if the two 's' planes form a 90 deg corner
tp@0:             planesare90degpairs = ismember( A(:,2),deg90planepairs(:,1) ).*ismember( A(:,3),deg90planepairs(:,2) ) + ismember( A(:,2),deg90planepairs(:,2) ).*ismember( A(:,3),deg90planepairs(:,1) );
tp@0:             ivsubset = find(planesare90degpairs);
tp@0:             ivdss = ivdss(ivsubset);
tp@0: 
tp@0:             if ~isempty(ivdss)
tp@0:                 addtocomb(ivdss,:) = [];
tp@0:                 addtooriginsfrom(ivdss,:) = [];
tp@0:                 addtoISESVISIBILITY(ivdss) = [];
tp@0:                 clear ivdss ivsubset planesare90degpairs                
tp@0:             end
tp@0:         end
tp@0:     end    
tp@0:     
tp@0:     nadds = size(addtocomb,1);
tp@0:     
tp@0:     if nadds > 0 & doconetrace == 1
tp@0: 
tp@0:         if difforder > 0
tp@0:             ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0:         else
tp@0:             ivss = uint32([1:nadds].');
tp@0:         end
tp@0:         
tp@0:         beamdirection1 = planemidpoints(addtocomb(ivss,ordernum-1),:) - ISCOORDS(addtooriginsfrom(ivss),1:3);
tp@0:         beamlength = sqrt( sum( beamdirection1.'.^2 ) ).';
tp@0:         beamdirection1  = beamdirection1./beamlength(:,ones(1,3));                
tp@0:         maxprojradius = maxdisttocorner(addtocomb(ivss,ordernum-1));
tp@0:         iv = find(beamlength<maxprojradius);
tp@0:         if ~isempty(iv)
tp@0:             beamlength(iv) = beamlength(iv)*0+ eps*1e6;    
tp@0:         end
tp@0:         beamangle1 = atan(maxprojradius./beamlength);
tp@0:                 
tp@0:         % Now we calculate the beam-widths to the second reflection
tp@0:         % plane, seen from the first-order image source. If this
tp@0:         % beamwidth is outside the previously calculated beam of
tp@0:         % reflection plane 1, then plane 2 can never be seen
tp@0:         % through plane 1.
tp@0: 
tp@0:         beamdirection2 = planemidpoints(addtocomb(ivss,ordernum),:) - ISCOORDS(addtooriginsfrom(ivss),1:3);
tp@0:         beamlength = sqrt( sum( beamdirection2.'.^2 ) ).';
tp@0:         beamdirection2  = beamdirection2./beamlength(:,ones(1,3));                
tp@0:         maxprojradius = maxdisttocorner(addtocomb(ivss,ordernum));
tp@0:         iv = find(beamlength<maxprojradius);
tp@0:         if ~isempty(iv)
tp@0:             beamlength(iv) = beamlength(iv)*0+ eps*10;    
tp@0:         end
tp@0:         beamangle2 = atan(maxprojradius./beamlength);
tp@0:        clear maxprojradius beamlength
tp@0:        
tp@0:        ivcutout = find( (beamangle1 + beamangle2 < acos(sum( (beamdirection1.*beamdirection2).' ).')) & (beamangle1 <1.56) & (beamangle2 < 1.56));
tp@0:        clear beamdirection1 beamdirection2 beamangle1 beamangle2
tp@0:        
tp@0:         if ~isempty(ivcutout)
tp@0: 			addtocomb(ivss(ivcutout),:) = [];
tp@0: 			addtooriginsfrom(ivss(ivcutout),:) = [];
tp@0:             addtoISESVISIBILITY(ivss(ivcutout)) = [];
tp@0:             nadds = size(addtocomb,1);
tp@0:         end
tp@0:                 
tp@0:     end
tp@0:     
tp@0: 	% Combinations with all-specular plus diffraction as the last one:
tp@0:     % check if the IS is behind both of the two planes that the reflecting edge connects to.
tp@0:     %
tp@0:     % Bug found 080711: The visibility test must be split into two, depending on the wedge angle!! 
tp@0:     % If the open wedge angle > pi, then it is correct to check if the IS behind both of the planes. 
tp@0:     % If the open wedge angle < pi, then it should be checked if the IS behind one of the planes!! 
tp@0: 
tp@0:     if difforder > 0
tp@0:         
tp@0:         ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0: 
tp@0:         if ~isempty(ivsd)
tp@0: 
tp@0:              % First we check the sssssd combinations where the open wedge angle > pi
tp@0:              % For those cases, we can exclude combinations where the IS
tp@0:              % is behind both planes
tp@0:                 
tp@0:              ivsd1 = ivsd( find(   closwedangvec( double(addtocomb(ivsd,ordernum))-nplanes ) < pi ));
tp@0: 
tp@0:              if ~isempty(ivsd1)
tp@0:                 imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1)    ,:).').';
tp@0:                 imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2)    ,:).').';
tp@0:                 GEOMACC = 1e-10;
tp@0:                 imsounotOK = find(imsoucheck1 < GEOMACC & imsoucheck2 < GEOMACC);
tp@0:                 clear imsoucheck1 imsoucheck2
tp@0:                 addtocomb(ivsd1(imsounotOK),:) = [];
tp@0:                 addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0:                 addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0:                 clear imsounotOK
tp@0:              end
tp@0:              
tp@0:              % Second we check the sssssd combinations where the open wedge angle < pi
tp@0:              % For those cases, we can exclude combinations where the IS
tp@0:              % is behind one of the two planes
tp@0:                 
tp@0:              ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0:              if ~isempty(ivsd)
tp@0:                  ivsd1 = ivsd( find(   closwedangvec( double(addtocomb(ivsd,ordernum))-nplanes ) > pi ));
tp@0: 
tp@0:                  if ~isempty(ivsd1)
tp@0:                     imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1)    ,:).').';
tp@0:                     imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2)     ,1),:)).',planenvecs(    planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2)    ,:).').';
tp@0:                     GEOMACC = 1e-10;
tp@0:                     imsounotOK = find(imsoucheck1 < GEOMACC | imsoucheck2 < GEOMACC);
tp@0:                     clear imsoucheck1 imsoucheck2
tp@0:                     addtocomb(ivsd1(imsounotOK),:) = [];
tp@0:                     addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0:                     addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0:                     clear imsounotOK
tp@0:                  end
tp@0:              end
tp@0:              
tp@0:             % For combinations including dsd, if the two diff. edges are
tp@0:             % aligned and the intermediate specular reflections are
tp@0:             % caused by planes that are perpendicular to the edges
tp@0:             % remove these combinations
tp@0:             
tp@0:             if difforder >= 2
tp@0:                 for kk = 1:ordernum-2
tp@0:                     ncombs = size(addtocomb,1); 
tp@0:                     matchvec = ones(ncombs,1);
tp@0:                     nprespecs = kk-1;
tp@0:                     nmidspecs = ordernum-1-kk;
tp@0:                     matchvec = matchvec.*(addtocomb(:,ordernum)>nplanes).*(addtocomb(:,kk)>nplanes);
tp@0:                     if nmidspecs == 1
tp@0:                         matchvec = matchvec.*(addtocomb(:,2)<=nplanes);
tp@0:                     elseif nmidspecs >= 2
tp@0:                         matchvec = matchvec.*prod( double(addtocomb(:,kk+1:ordernum-1)<=nplanes).' ).';
tp@0:                     end
tp@0:                     ivdsd = uint32(find( matchvec ));
tp@0:                     if ~isempty(ivdsd)
tp@0:                         edge1 = double(addtocomb(ivdsd,kk)) - nplanes;
tp@0:                         edge2 = double(addtocomb(ivdsd,ordernum)) - nplanes;
tp@0:                         midspecs = double(addtocomb(ivdsd,kk+1:ordernum-1));
tp@0:                         
tp@0:                         lookupindmat1 = (edge1-1)*nedges + edge2;
tp@0:                         matrixcolnumbers = (edge1-1)*nplanes;
tp@0:                         lookupindmat2 = matrixcolnumbers(:,ones(1,nmidspecs)) + midspecs;                    
tp@0:                         
tp@0:                         if nmidspecs == 1
tp@0:                             ivinvalid = find(edgealignedwithedge(lookupindmat1) & edgeperptoplane(lookupindmat2));
tp@0:                         else
tp@0:                             ivinvalid = find(edgealignedwithedge(lookupindmat1) & prod(edgeperptoplane(lookupindmat2).').');
tp@0:                         end
tp@0:             			addtocomb(ivdsd(ivinvalid),:) = [];
tp@0: 				        addtooriginsfrom(ivdsd(ivinvalid),:) = [];
tp@0:                         addtoISESVISIBILITY(ivdsd(ivinvalid)) = [];
tp@0:                     end
tp@0:                  end
tp@0:             end            
tp@0:             
tp@0:             % We need to find the valid combinations again after we have removed a
tp@0:             % number of combinations
tp@0: 	
tp@0:             ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0: 	
tp@0:             %----------------------------------------------------------------------
tp@0:             % Combinations of specular-diffraction that are not visible at
tp@0:             % all should be removed and then they are not propagated either.
tp@0: 		
tp@0:             if nhighval < 255
tp@0:                 possibleedges = uint8(double(addtocomb(ivsd,ordernum)) - nplanes);
tp@0:                 possiblecombs = uint8(double(addtocomb(ivsd,1:ordernum-1)));
tp@0:             else
tp@0:                 possibleedges = uint16(double(addtocomb(ivsd,ordernum)) - nplanes);
tp@0:                 possiblecombs = uint16(double(addtocomb(ivsd,1:ordernum-1)));
tp@0:             end
tp@0:             reftoISCOORDS = addtooriginsfrom(ivsd);
tp@0: 		
tp@0:             % Expand to take the edge subdivisions into account
tp@0:             
tp@0:             nposs = length(ivsd);
tp@0:             nposs = nposs*nedgesubs;        % We treat all the little edge subdivisions as separate edges
tp@0: 
tp@0:             expandedposscombs = reshape(repmat(possiblecombs.',[nedgesubs,1]),ordernum-1,nposs).';
tp@0:             clear possiblecombs                        
tp@0:             expandedreftoISCOORDS = reftoISCOORDS(:,onesvec1);
tp@0:             expandedreftoISCOORDS = reshape(expandedreftoISCOORDS.',nposs,1);
tp@0:             expandedpossedges = possibleedges(:,onesvec1);
tp@0:             expandedpossedges = reshape(expandedpossedges.',nposs,1);
tp@0:             expandedivsd = ivsd(:,onesvec1);
tp@0:             expandedivsd = reshape(expandedivsd.',nposs,1);
tp@0:             
tp@0:             if SHOWTEXT >= 3
tp@0: 				disp(['         ',int2str(nposs),' IS+edge segments found initially '])
tp@0: 			end
tp@0: 	
tp@0:             % Go through, iteratively, and check if the path from S to edge passes
tp@0:             % through all reflection planes along the way
tp@0:             
tp@0:             for jj = ordernum-1:-1:1
tp@0: 	
tp@0:                 if nposs > 0
tp@0: 		
tp@0:                     if jj == ordernum-1
tp@0:                         fromcoords = ISCOORDS(expandedreftoISCOORDS,:);
tp@0:                         [tocoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0:                         clear possibleedges
tp@0:                     else
tp@0:                         eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])    
tp@0:                         ivlist = ORIGINSFROM(expandedreftoISCOORDS);
tp@0:                         for kk = jj:ordernum-3
tp@0:                             ivlist = ORIGINSFROM(ivlist);                            
tp@0:                         end
tp@0:                         fromcoords = ISCOORDS(ivlist,:);
tp@0:                                            
tp@0:                     end
tp@0: 
tp@0:                     [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,jj),4),planenvecs(expandedposscombs(:,jj),:),minvals(expandedposscombs(:,jj),:),...
tp@0: 					    maxvals(expandedposscombs(:,jj),:),planecorners(expandedposscombs(:,jj),:),corners,ncornersperplanevec(expandedposscombs(:,jj)));
tp@0:                     if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0:                         disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0:                         disp('         handled correctly yet.')
tp@0:                     end
tp@0:                     eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0: 			
tp@0:                     expandedivsd          = expandedivsd(hitplanes);
tp@0:                     expandedposscombs     = expandedposscombs(hitplanes,:);
tp@0:                     expandedreftoISCOORDS = expandedreftoISCOORDS(hitplanes);
tp@0:                     expandedpossedges     = expandedpossedges(hitplanes);
tp@0:                     edgeweightlist        = edgeweightlist(hitplanes);
tp@0:                     toedgecoords          = tocoords(hitplanes,:);
tp@0:                     if jj < ordernum-1
tp@0:                         for kk = jj+1:ordernum-1
tp@0:                            eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);'])
tp@0:                         end
tp@0:                     end
tp@0: 		
tp@0:                     nposs = length(expandedivsd);
tp@0:                     
tp@0:                 end
tp@0:                 if SHOWTEXT >= 3
tp@0: 				    disp(['         ',int2str(nposs),' IS+edge segments survived the visibility test in refl plane ',JJ(jj,1:JJnumbofchars(jj))])
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0: 	
tp@0:             if obstructtestneeded & nposs > 0
tp@0:                 for jj = 1:ordernum
tp@0:                     if nposs > 0
tp@0:                         
tp@0:                         if jj==1
tp@0:                             fromcoords = S;    
tp@0:                             startplanes = [];    
tp@0:                         else
tp@0:                             startplanes = expandedposscombs(:,jj-1);
tp@0:                             eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0:                         end
tp@0:                         if jj == ordernum
tp@0:                             tocoords = toedgecoords;
tp@0:                             endplanes = [planesatedge(expandedpossedges,1) planesatedge(expandedpossedges,2)];
tp@0:                         else
tp@0:                             eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])    
tp@0:                             endplanes = expandedposscombs(:,jj);    
tp@0:                         end
tp@0:              
tp@0:                         [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0:                             planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0: 		
tp@0:                         if nobstructions > 0
tp@0:                             expandedivsd          = expandedivsd(nonobstructedpaths);
tp@0:                             expandedposscombs     = expandedposscombs(nonobstructedpaths,:);
tp@0:                             expandedreftoISCOORDS = expandedreftoISCOORDS(nonobstructedpaths);
tp@0:                             expandedpossedges     = expandedpossedges(nonobstructedpaths);
tp@0:                             edgeweightlist        = edgeweightlist(nonobstructedpaths);
tp@0:                             toedgecoords          = tocoords(nonobstructedpaths,:);
tp@0:                             nposs = length(expandedivsd);
tp@0:                             for kk = 1:ordernum-1
tp@0:                                 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])    
tp@0:                             end
tp@0:                             
tp@0:                         end
tp@0:                     end
tp@0:                 end
tp@0:             end        
tp@0:         
tp@0:             if SHOWTEXT >= 3
tp@0: 		        disp(['         ',int2str(nposs),' IS+edge segments survived the obstruction test'])
tp@0:             end
tp@0: 
tp@0:             % There are repetitions of each combination since each edge segment has
tp@0:             % been treated separately. Pack them together now
tp@0: 	
tp@0:             test = [expandedposscombs expandedpossedges];
tp@0:             ncombs = length(expandedpossedges);
tp@0:             dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
tp@0:             ivremove = find(dtest==1);
tp@0:             
tp@0:             while ~isempty(ivremove)
tp@0:            
tp@0:                 edgeweightlist(ivremove+1) = double(edgeweightlist(ivremove+1)) + double(edgeweightlist(ivremove));
tp@0:                 edgeweightlist(ivremove) = [];
tp@0:                 expandedpossedges(ivremove) = [];
tp@0:                 expandedposscombs(ivremove,:) = [];
tp@0:                 expandedivsd(ivremove) = [];
tp@0:                 nposs = length(expandedivsd);   
tp@0:                 
tp@0:                 test = [expandedposscombs expandedpossedges];
tp@0:                 ncombs = length(expandedpossedges);
tp@0:                 dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
tp@0:                 ivremove = find(dtest==1);
tp@0: 	
tp@0:             end
tp@0:             	       
tp@0:             if SHOWTEXT >= 3
tp@0: 				disp(['         ',int2str(nposs),' IS+edge segments after edge segments are packed together'])
tp@0:             end
tp@0: 	
tp@0:             combstoremove = setdiff(ivsd,expandedivsd);
tp@0:             addtocomb(combstoremove,:) = [];
tp@0:             addtooriginsfrom(combstoremove) = [];  
tp@0:             addtoISESVISIBILITY(combstoremove) = [];
tp@0:             nadds = length(addtooriginsfrom);
tp@0:         end        
tp@0:     end   
tp@0: 
tp@0:     POTENTIALISES = [POTENTIALISES;[addtocomb zeros(nadds,specorder-ordernum)]];
tp@0:     ORIGINSFROM = [ORIGINSFROM;addtooriginsfrom];
tp@0: 
tp@0:     ivtemp = find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes));
tp@0:     if difforder > 0
tp@0:         if ~isempty(ivtemp)
tp@0:             addtoISESVISIBILITY(ivtemp) = edgeweightlist;
tp@0:         end
tp@0:     end
tp@0:     ISESVISIBILITY = [ISESVISIBILITY;addtoISESVISIBILITY];
tp@0:     
tp@0:     endindices(ordernum) = length(ORIGINSFROM);
tp@0: 
tp@0:     % Compute the IS coords
tp@0: 
tp@0:     ISCOORDStoadd = zeros(nadds,3);
tp@0:     if difforder > 0
tp@0:         ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0:     else
tp@0:         ivss = uint32([1:nadds].');    
tp@0:     end
tp@0:     soutoreflect = ISCOORDS(ORIGINSFROM(double(ivss)+endindices(ordernum-1)),1:3);
tp@0:     ISCOORDStoadd(ivss,:) = EDB1findis(soutoreflect,POTENTIALISES(double(ivss)+endindices(ordernum-1),ordernum),planeeqs,size(soutoreflect,1),onesvec3);
tp@0:     clear soutoreflect
tp@0:     ISCOORDS = [ISCOORDS;ISCOORDStoadd];
tp@0:     clear ISCOORDStoadd
tp@0:     
tp@0: end
tp@0: 
tp@0: ntot = size(POTENTIALISES,1);
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Add the direct sound to the start of the list
tp@0: 
tp@0: POTENTIALISES = [zeros(1,specorder);POTENTIALISES];
tp@0: ORIGINSFROM = [0;ORIGINSFROM];
tp@0: startindices = startindices+1;
tp@0: endindices = endindices+1;
tp@0: ORIGINSFROM = uint32(double(ORIGINSFROM)+1);
tp@0: 
tp@0: ISCOORDS = [S;ISCOORDS];
tp@0: ISESVISIBILITY = [0;ISESVISIBILITY];
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Trim the list if there are zeros-only in the last column(s)
tp@0: 
tp@0: [n1,n2] = size(POTENTIALISES);
tp@0: if n2 > 1
tp@0:     columnstokeep = find(sum(POTENTIALISES)~=0);
tp@0:     POTENTIALISES = POTENTIALISES(:,columnstokeep);    
tp@0:     [n1,n2new] = size(POTENTIALISES);    
tp@0:     if n2new < n2
tp@0:         specorderold = specorder;
tp@0:         specorder = n2new;
tp@0:     end
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Create index lists
tp@0: 
tp@0: % First, the total reflection order
tp@0: 
tp@0: [n1,n2] = size(POTENTIALISES);
tp@0: 
tp@0: lengthNspecmatrix = [];
tp@0: lengthNdiffmatrix = [];
tp@0: singlediffcol = [];
tp@0: startindicessinglediff = [];
tp@0: endindicessinglediff = [];
tp@0: 
tp@0: if n1 > 0 & n2 > 0
tp@0: if n2 > 1
tp@0:     REFLORDER = uint8(sum(POTENTIALISES.'>0).');    
tp@0: else
tp@0:     REFLORDER = uint8(POTENTIALISES>0);
tp@0: end
tp@0: 
tp@0: if difforder > 0
tp@0:     if specorder > 1
tp@0:         ivss = uint32(find(prod( double(POTENTIALISES<=nplanes).' ).'));
tp@0:     else
tp@0:         ivss = uint32(find( POTENTIALISES<=nplanes ));    
tp@0:     end
tp@0: 
tp@0:     ivother = uint32([1:length(ORIGINSFROM)].');
tp@0:     ivother(ivss) = [];
tp@0:     ivss(1) = [];
tp@0:     IVNDIFFMATRIX = uint32(zeros(2,specorder));
tp@0:     lengthNdiffmatrix = zeros(1,specorder);
tp@0:     nrowsdiff = 2;
tp@0:     IVNSPECMATRIX = uint32(zeros(2,specorder));
tp@0:     lengthNspecmatrix = zeros(1,specorder);
tp@0:     nrowsspec = 2;
tp@0:     
tp@0:     for ii = 1:specorder
tp@0:         if specorder > 1
tp@0:             ivreftoNdiff = uint32(find(sum( double(POTENTIALISES(ivother,:)> nplanes).' ).'==ii));
tp@0:         else
tp@0:             ivreftoNdiff = uint32(find(( double(POTENTIALISES(ivother,:)> nplanes).' ).'==ii));
tp@0:         end
tp@0:         ivreftoNspec = uint32(find(REFLORDER(ivss)==ii));
tp@0: 
tp@0:         ivNdiff = ivother(ivreftoNdiff);
tp@0:         if ~isempty(ivNdiff)
tp@0:             lengthNdiffmatrix(ii) = length(ivNdiff);
tp@0:             if lengthNdiffmatrix(ii) > nrowsdiff
tp@0:                IVNDIFFMATRIX = [IVNDIFFMATRIX;uint32(zeros(lengthNdiffmatrix(ii)-nrowsdiff,specorder))];
tp@0:                nrowsdiff = lengthNdiffmatrix(ii);
tp@0:             end
tp@0:             IVNDIFFMATRIX(1: lengthNdiffmatrix(ii),ii) =  ivNdiff;     
tp@0:         end
tp@0:         ivother(ivreftoNdiff) = [];
tp@0:         
tp@0:         ivNspec = ivss(ivreftoNspec);
tp@0:         if ~isempty(ivNspec)
tp@0:             lengthNspecmatrix(ii) = length(ivNspec);
tp@0:             if lengthNspecmatrix(ii) > nrowsspec
tp@0:                 IVNSPECMATRIX = [IVNSPECMATRIX;uint32(zeros(lengthNspecmatrix(ii)-nrowsspec,specorder))];
tp@0:                 nrowsspec = lengthNspecmatrix(ii);
tp@0:             end
tp@0:             IVNSPECMATRIX(1: lengthNspecmatrix(ii),ii) =  ivNspec;      
tp@0:         end        
tp@0:     end
tp@0:      
tp@0:     % Determine which column the single diffraction is in
tp@0:     
tp@0:     test = POTENTIALISES(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1),:)>nplanes;
tp@0:        
tp@0:     colweight = [1:specorder];
tp@0:     nsingles = length(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1));
tp@0:     if specorder > 1
tp@0:         singlediffcol = uint8(sum( (test.*colweight(ones(nsingles,1),:)).').');
tp@0:     else
tp@0:         singlediffcol = uint8(ones(nsingles,1));        
tp@0:     end
tp@0:         
tp@0:     startindicessinglediff = zeros(specorder,1);
tp@0:     endindicessinglediff = zeros(specorder,1);
tp@0:     startindicessinglediff(1) = 1;
tp@0:     for ii = 1:specorder-1
tp@0:         iv = find(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1)<=endindices(ii));
tp@0:         if ~isempty(iv)
tp@0:             endindicessinglediff(ii) = iv(length(iv));
tp@0:         else
tp@0:             endindicessinglediff(ii) = 0;
tp@0:         end
tp@0:         startindicessinglediff(ii+1) = endindicessinglediff(ii)+1;
tp@0:     end
tp@0:     endindicessinglediff(specorder) = lengthNdiffmatrix(1);
tp@0: else
tp@0:     ivss = uint32([2:endindices(length(endindices))].');
tp@0:     IVNSPECMATRIX = uint32(zeros(2,specorder));
tp@0:     lengthNspecmatrix = zeros(1,specorder);
tp@0:     nrowsspec = 2;
tp@0: 
tp@0:     for ii = 1:specorder
tp@0:         ivreftoNspec = uint32(find(REFLORDER(ivss)==ii));
tp@0: 
tp@0:         ivNspec = ivss(ivreftoNspec);
tp@0:         lengthNspecmatrix(ii) = length(ivNspec);
tp@0:         if lengthNspecmatrix(ii) > nrowsspec
tp@0:            IVNSPECMATRIX = [IVNSPECMATRIX;uint32(zeros(lengthNspecmatrix(ii)-nrowsspec,specorder))];
tp@0:            nrowsspec = lengthNspecmatrix(ii);
tp@0:         end
tp@0:         IVNSPECMATRIX(1: lengthNspecmatrix(ii),ii) =  ivNspec;      
tp@0:     
tp@0:     end
tp@0:     
tp@0:     IVNDIFFMATRIX = [];
tp@0:     lengthNdiffmatrix = [];
tp@0: 
tp@0:     singlediffcol = [];
tp@0:     startindicessinglediff = [];
tp@0:     endindicessinglediff = [];
tp@0: end
tp@0: 
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Create a pointer list so that it is possible to find a combination
tp@0: % that ends with edge-plane2, given an edge-plane1-plane2 combination
tp@0: %
tp@0: % NB!! The list points to the original POTENTIALISES (and related lists)
tp@0: 
tp@0: if difforder > 0 & specorder > 1
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0:         disp(['       Building a pointer list for image receiver combinations'])    
tp@0:     end
tp@0: 
tp@0:     % First select only those that have a single diff combo and that
tp@0:     % start with the diffraction
tp@0:     
tp@0:     ndecimaldivider = (nplanes+2);
tp@0:     PointertoIRcombs = sparse(zeros(ndecimaldivider^(specorder-1),1));
tp@0:     IRoriginsfrom = uint32(zeros(size(ORIGINSFROM)));
tp@0:     
tp@0: 	for ii = 1:specorder-1
tp@0:         if SHOWTEXT >= 3
tp@0:             disp(['          order ',int2str(ii)])    
tp@0:         end
tp@0:                 
tp@0:         if ii == 1
tp@0:             ivrange = uint32([startindicessinglediff(2):endindicessinglediff(2)]);
tp@0:             masterivlistselect = IVNDIFFMATRIX(ivrange,1);
tp@0:             masterivlistselect = masterivlistselect(find(singlediffcol(ivrange)==1));
tp@0: 
tp@0:             % IRoriginsfrom should be given the value zero for these
tp@0:             % first-order specular combos (after one diff) so we don't
tp@0:             % bother assigning any value
tp@0:             
tp@0:             [B,IA,JA] = unique(POTENTIALISES(masterivlistselect,2));
tp@0:             PointertoIRcombs(B) = masterivlistselect(IA);
tp@0:             
tp@0:             % Now we check if there any active wall numbers that don't
tp@0:             % occur in an edge-spec combination. For those combinations
tp@0:             % we can point to a specular combos instead, that ends with the
tp@0:             % same plane number.
tp@0:             
tp@0:             wallist = POTENTIALISES(IVNSPECMATRIX(1:lengthNspecmatrix(1),1),1);
tp@0:             ivnotincludedyet = find(~ismember(wallist,B));
tp@0:             if ~isempty(ivnotincludedyet)
tp@0:                 PointertoIRcombs(wallist(ivnotincludedyet)) = IVNSPECMATRIX(ivnotincludedyet,1);
tp@0:             end
tp@0:             
tp@0:         elseif ii >= 2
tp@0:             ivrange = uint32([startindicessinglediff(ii+1):endindicessinglediff(ii+1)]);
tp@0:             masterivlistselect = IVNDIFFMATRIX(ivrange,1);
tp@0:             masterivlistselect = masterivlistselect(find(singlediffcol(ivrange)==1));
tp@0:  
tp@0:             ivlist = 0;
tp@0:             for jj = 3:ii+1
tp@0:                 ivlist = ivlist + double(POTENTIALISES(masterivlistselect,jj))*ndecimaldivider^(ii+1-jj);
tp@0:             end
tp@0:             
tp@0:             IRoriginsfrom(masterivlistselect) = uint32(full(PointertoIRcombs(ivlist)));
tp@0: 
tp@0:             ivlist = ivlist + double(POTENTIALISES(masterivlistselect,2))*ndecimaldivider^(ii-1);
tp@0:             A = uint32(ivlist);
tp@0:             [B,IA,JA] = unique(A);
tp@0:             PointertoIRcombs(B) = masterivlistselect(IA);
tp@0: 
tp@0:             % Now we check if there any active wall numbers that don't
tp@0:             % occur as spec1 in an edge-spec1-spec2-spec3 combination.
tp@0: 
tp@0:             wallist = 0;
tp@0:             for jj = 1:ii
tp@0:                 wallist = wallist + double(POTENTIALISES(IVNSPECMATRIX(1:lengthNspecmatrix(ii),ii),jj))*ndecimaldivider^(ii-jj);
tp@0:             end
tp@0: 
tp@0:             ivnotincludedyet = find(~ismember(wallist,B));
tp@0: 
tp@0:             if ~isempty(ivnotincludedyet)
tp@0:                 PointertoIRcombs(wallist(ivnotincludedyet)) = IVNSPECMATRIX(ivnotincludedyet,ii);
tp@0:             end
tp@0:             
tp@0:         end
tp@0: 	
tp@0: 	end
tp@0:         
tp@0: else
tp@0:     ndecimaldivider = 0;
tp@0:     PointertoIRcombs = [];
tp@0:     IRoriginsfrom = [];
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Save the output data
tp@0: 
tp@0: maxval = max(IRoriginsfrom);
tp@0: if maxval < 256
tp@0:     IRoriginsfrom = uint8(IRoriginsfrom);
tp@0: elseif maxval < 65536
tp@0:     IRoriginsfrom = uint16(IRoriginsfrom);    
tp@0: end
tp@0: maxval = max(ORIGINSFROM);
tp@0: if maxval < 256
tp@0:     ORIGINSFROM = uint8(ORIGINSFROM);
tp@0: elseif maxval < 65536
tp@0:     ORIGINSFROM = uint16(ORIGINSFROM);    
tp@0: end
tp@0: 
tp@0: if SUPPRESSFILES == 0
tp@0:     % The complete variable set could be saved if one wanted to.
tp@0:     Varlist = [' POTENTIALISES ORIGINSFROM ISCOORDS  ISESVISIBILITY IVNSPECMATRIX lengthNspecmatrix IVNDIFFMATRIX lengthNdiffmatrix '];
tp@0:     Varlist = [Varlist,' singlediffcol REFLORDER startindicessinglediff endindicessinglediff ndecimaldivider PointertoIRcombs IRoriginsfrom'];		
tp@0: 
tp@0:     % Varlist = [' lengthNspecmatrix lengthNdiffmatrix '];
tp@0:     % Varlist = [Varlist,' singlediffcol startindicessinglediff endindicessinglediff ndecimaldivider PointertoIRcombs IRoriginsfrom'];		
tp@0: 
tp@0:     eval(['save ',ISEStreefile,Varlist])
tp@0: else
tp@0:    ISEStreefile = struct('POTENTIALISES',POTENTIALISES,'ORIGINSFROM',ORIGINSFROM,'ISCOORDS',ISCOORDS,'ISESVISIBILITY',ISESVISIBILITY,'IVNSPECMATRIX',IVNSPECMATRIX,...
tp@0:        'lengthNspecmatrix',lengthNspecmatrix,'IVNDIFFMATRIX',IVNDIFFMATRIX,'lengthNdiffmatrix',lengthNdiffmatrix,'singlediffcol',singlediffcol,...
tp@0:        'REFLORDER',REFLORDER,'startindicessinglediff',startindicessinglediff,'endindicessinglediff',endindicessinglediff,'ndecimaldivider',ndecimaldivider,...
tp@0:        'PointertoIRcombs',PointertoIRcombs,'IRoriginsfrom',IRoriginsfrom);    
tp@0: end
tp@0:     
tp@0:     
tp@0:     
tp@0: