tp@0: function [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,listguide,listoforders,...
tp@0:         bigedgeweightlist] = EDB1diffISESx(eddatafile,S,R,...
tp@0:         ivsinglediff,singlediffcol,startindicessinglediff,endindicessinglediff,...
tp@0:         specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,...
tp@0:         PointertoIRcombs,IRoriginsfrom)
tp@0: % EDB1diffISESx - Gives list of paths that includes a 1st-order diff. combination.
tp@0: % Gives the list of possible first-order diffraction paths, possibly with specular
tp@0: % reflections before and after.
tp@0: %
tp@0: % Input parameters:
tp@0: %   eddatafile      The name of the file that contains all edge related data.
tp@0: %                   This file will be loaded.
tp@0: %   S,R,ivsinglediff,...
tp@0: %    singlediffcol,startindicessinglediff,endindicessinglediff,specorder,visplanesfromR,...
tp@0: %   vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,PointertoIRcombs,IRoriginsfrom
tp@0: %                   Data that should have been passed from the srdatafile
tp@0: %                   (S,R,visplanesfromR,vispartedgesfromS,vispartedgesfromR
tp@0: %                   )
tp@0: %                   from the ISEStreefile
tp@0: %                   (ivsinglediff
tp@0: %                   singlediffcol,startindicessinglediff,endindicessinglediff
tp@0: %                   ndecimaldivider,PointertoIRcombs,IRoriginsfrom)
tp@0: %                   and from the setupfile
tp@0: %                   (specorder,nedgesubs)
tp@0: %   POTENTIALISES (global)
tp@0: %   ISCOORDS (global)
tp@0: %   ORIGINSFROM (global)
tp@0: %   ISESVISIBILITY (global)   
tp@0: %   REFLORDER (global)
tp@0: %
tp@0: % Global parameters:
tp@0: %   SHOWTEXT,JJ,JJnumbofchars     See EDB1mainISES
tp@0: %   POTENTIALISES,ISCOORDS
tp@0: %
tp@0: % Output parameters:
tp@0: %   edgedifflist        List [ncombs,1] of the edge number involved in each
tp@0: %                       spec-diff-spec combination.
tp@0: %   startandendpoints   Matrix [ncombs,2] of the relative start and end
tp@0: %                       points of each edge. The values, [0,1], indicate
tp@0: %                       which part of the edge that is visible.
tp@0: %   prespeclist         Matrix [ncombs,specorder-1] of the specular
tp@0: %                       reflections that precede every diffraction.
tp@0: %   postspeclist        Matrix [ncombs,specorder-1] of the specular
tp@0: %                       reflections that follow every diffraction.
tp@0: %   validIScoords       Matrix [ncombs,3] of the image source for each
tp@0: %                       multiple-spec that precedes the diffraction. If
tp@0: %                       there is no spec refl before the diffraction, the
tp@0: %                       value [0 0 0] is given.
tp@0: %   validIRcoords       Matrix [ncombs,3] of the image receiver for each
tp@0: %                       multiple-spec that follows the diffraction. If
tp@0: %                       there is no spec refl after the diffraction, the
tp@0: %                       value [0 0 0] is given.
tp@0: %   listguide           Matrix [nuniquecombs,3] which for each row gives
tp@0: %                       1. The number of examples in edgefdifflist etc that
tp@0: %                          are the same type of spec-diff-spec comb.
tp@0: %                       2. The first row number and 3. The last row number.
tp@0: %   listoforders        Matrix [nuniquecombs,2] which for each row gives
tp@0: %                       1. The reflection order for the spec-diff-spec comb
tp@0: %                          in the same row in listguide.
tp@0: %                       2. The order of the diffraction in the
tp@0: %                          spec-diff-spec comb.
tp@0: %   bigedgeweightlist   List [ncombs,1] of the visibility of each edge
tp@0: %                       expressed as a number 0 to 2^nedgesubs-1.
tp@0: %
tp@0: % Uses functions  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) 20061118
tp@0: %
tp@0: % [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,listguide,listoforders,...
tp@0: %   bigedgeweightlist] = EDB1diffISESx(eddatafile,S,R,...
tp@0: %   ivsinglediff,singlediffcol,startindicessinglediff,endindicessinglediff,...
tp@0: %   specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,...
tp@0: %   PointertoIRcombs,IRoriginsfrom)
tp@0: 
tp@0: global SHOWTEXT JJ JJnumbofchars
tp@0: global POTENTIALISES ISCOORDS ORIGINSFROM REFLORDER ISESVISIBILITY
tp@0: 
tp@0: eval(['load ',eddatafile])
tp@0: clear edgeseesplane cornerinfrontofplane
tp@0: 
tp@0: [nedges,slask] = size(planesatedge);
tp@0: [nplanes,slask] = size(planecorners);
tp@0: 
tp@0: edgedifflist = [];
tp@0: prespeclist = [];
tp@0: postspeclist = [];
tp@0: startandendpoints = [];
tp@0: bigedgeweightlist = [];
tp@0: validIScoords = [];
tp@0: validIRcoords = [];
tp@0: 
tp@0: [n1,n2] = size(POTENTIALISES);
tp@0: if n2 < specorder
tp@0:     specorder = n2;    
tp@0: end
tp@0: 
tp@0: maxvisibilityvalue = 2^nedgesubs-1;
tp@0: zerosvec1 = zeros(1,max([specorder-1 1]));
tp@0: zerosvec2 = zeros(1,3);
tp@0: listguide = zeros(specorder*2-1,3);
tp@0: listoforders = zeros(specorder*2-1,2);
tp@0: 
tp@0: obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0: onesvec = ones(1,nedgesubs);
tp@0: onesvec3 = ones(1,3);
tp@0: 
tp@0: %   ###########################################
tp@0: %   #                                         #
tp@0: %   #         S - spec-spec- edge - R cases   #
tp@0: %   #                                         #
tp@0: %   #         Prespec cases                   #
tp@0: %   #                                         #
tp@0: %   ###########################################
tp@0: %
tp@0: % Possible edges for S-spec-spec-E-R are seen (at least partly) by the receiver.
tp@0: %
tp@0: % The visibility doesn't need to be checked since this was done in the
tp@0: % ISEStree.
tp@0: 
tp@0: % The vector masterivlist will always refer to the original data vector
tp@0: % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc
tp@0: %
tp@0: % First we pick out those indices where there was a single diffraction, but
tp@0: % skip those with only diffraction (because we dealt with them already).
tp@0: % Also, select only those where the diffraction is the last in the sequence
tp@0: % of reflections.
tp@0: 
tp@0: for ii = 1:specorder,       % NB!!! ii = 1 corresponds to zero spec. refl before the diff.
tp@0:                             %       ii = 2 corresponds to one spec refl.
tp@0:                             %       before the diff etc
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0:         if ii > 1
tp@0:             disp(['      Checking for ',JJ(ii-1,1:JJnumbofchars(ii-1)),' spec refl before the edge diff'])    
tp@0:         else
tp@0:             disp(['      Checking for 0 spec refl before the edge diff'])    
tp@0:         end
tp@0:     end
tp@0: 
tp@0:     iv = uint32(startindicessinglediff(ii):endindicessinglediff(ii));
tp@0:         
tp@0:     if ~isempty(iv)
tp@0:                 
tp@0:         ivkeep = find(singlediffcol(iv)==ii);
tp@0:         iv = iv(ivkeep);
tp@0:         masterivlist = ivsinglediff(iv);
tp@0:         possibleedges = double(POTENTIALISES(masterivlist,ii)) - nplanes;
tp@0: 
tp@0:         % Keep only combinations for which the receiver can see the edge
tp@0:             
tp@0:         ivnotvisiblefromr = find(vispartedgesfromR(possibleedges)==0);
tp@0:         if ~isempty(ivnotvisiblefromr)
tp@0:             masterivlist(ivnotvisiblefromr) = [];
tp@0:             possibleedges(ivnotvisiblefromr) = [];
tp@0:         end        
tp@0:         % Pick out the pre-specs
tp@0: 
tp@0:         nposs = length(masterivlist);
tp@0: 
tp@0:         if nposs > 0
tp@0:         
tp@0:             if ii > 1
tp@0:                 possiblecombs = POTENTIALISES(masterivlist,1:ii-1);
tp@0:                 reftoIScoords = ORIGINSFROM(masterivlist);
tp@0:                 edgeweightlist = bitand((ISESVISIBILITY(masterivlist)),(vispartedgesfromR(possibleedges)));
tp@0:                 prespeclist = [prespeclist;[possiblecombs zeros(nposs,specorder-ii)]];
tp@0:                 % NB! It is correct below that the indices for the IScoords should be
tp@0:                 % ORIGINSFROM(masterivlist), rather than masterivlist.
tp@0:                 % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0:                 % spec-spec-...-diff combinations and then
tp@0:                 % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0:                 % ends with a diff has no image source. 
tp@0:                  validIScoords = [validIScoords;ISCOORDS(ORIGINSFROM(masterivlist),:)];
tp@0:             else            
tp@0:                 edgeweightlist = bitand((vispartedgesfromS(possibleedges)),(vispartedgesfromR(possibleedges)));
tp@0:                 prespeclist = [prespeclist;[zeros(nposs,max([specorder-1 1]))]];       
tp@0:                 % For the case of no spec refl before the diffraction, we
tp@0:                 % let the IS get the S coordinates.
tp@0:                 validIScoords = [validIScoords;S(ones(nposs,1),:)];
tp@0:             end
tp@0:             
tp@0:             if SHOWTEXT >= 3
tp@0:          		disp(['         ',int2str(nposs),' IS+edges valid'])
tp@0:             end
tp@0: 		
tp@0:             edgedifflist = [edgedifflist;possibleedges];
tp@0:             postspeclist = [postspeclist;zerosvec1(ones(nposs,1),:)];
tp@0:             bigedgeweightlist = [bigedgeweightlist;edgeweightlist];   
tp@0:             validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0:         
tp@0:         end
tp@0:         
tp@0:     end
tp@0: end
tp@0: 
tp@0: iv = uint32(find(bigedgeweightlist==0));
tp@0: if ~isempty(iv)
tp@0:     edgedifflist(iv) = [];
tp@0:     bigedgeweightlist(iv) = [];
tp@0:     prespeclist(iv,:) = [];
tp@0:     postspeclist(iv,:) = [];
tp@0:     validIScoords(iv,:) = [];
tp@0:     validIRcoords(iv,:) = [];
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: 
tp@0: %   #############################################
tp@0: %   #                                           #
tp@0: %   #         Build an IR tree                 #
tp@0: %   #                                           #
tp@0: %   #############################################
tp@0: %
tp@0: % For all cases with a specular reflection after an edge diffraction
tp@0: % we build an IR tree, analogous to the IS tree
tp@0: 
tp@0: if specorder > 1
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0:             disp(['      For the edge+spec combs, an IR tree is built'])    
tp@0: 	end
tp@0: 
tp@0:     % Select all combinations with a single diffraction, anywhere except in
tp@0:     % the last column, and store the useful data for these (which column
tp@0:     % is the diffraction, what reflection order, how many specular
tp@0:     % reflections?)
tp@0:    
tp@0:     ivselect = uint32(find( singlediffcol~=REFLORDER(ivsinglediff) ));
tp@0:     masterivlistorig = ivsinglediff(ivselect);
tp@0:     singlediffcol_select = singlediffcol(ivselect);
tp@0:     reflorder_select = REFLORDER(ivsinglediff(ivselect));
tp@0:     clear ivsinglediff
tp@0:     nspecular_select = uint32(double(reflorder_select) - double(singlediffcol_select));
tp@0:     nspecular_select_not_empty = 1;
tp@0:     
tp@0:     % Now we remove all those with a last reflection plane which can not be
tp@0:     % seen by the receiver
tp@0:     
tp@0:     if nplanes + nedges < 65536
tp@0:         lastreflplane = uint16(zeros(length(ivselect),1));
tp@0:     else
tp@0:         lastreflplane = uint32(zeros(length(ivselect),1));        
tp@0:     end
tp@0:     clear ivselect
tp@0:     for ii = 2:specorder
tp@0:         iv = find(reflorder_select==ii);
tp@0:         lastreflplane(iv) = POTENTIALISES(masterivlistorig(iv),ii);
tp@0:     end
tp@0:     ivremove = uint32(find(visplanesfromR(lastreflplane)~=2 & visplanesfromR(lastreflplane)~=4 & visplanesfromR(lastreflplane)~=5));
tp@0:     if ~isempty(ivremove)
tp@0:         masterivlistorig(ivremove) = [];
tp@0:         singlediffcol_select(ivremove) = [];
tp@0:         reflorder_select(ivremove) = [];
tp@0:         nspecular_select(ivremove) = [];
tp@0:         lastreflplane(ivremove) = [];
tp@0:         clear ivremove
tp@0:     end  
tp@0:     	
tp@0:     % Start by calculating all the first-order IR coordinates for
tp@0:     % all the planes that are visible from the receiver
tp@0: 	IRcoordsstart = zeros(nplanes,3);
tp@0: 	iv = uint32(find(visplanesfromR==2 | visplanesfromR==4 | visplanesfromR==5));
tp@0: 	IRcoordsstart(iv,:) = EDB1findis(R,iv,planeeqs,1,onesvec3);
tp@0: 		
tp@0: 	bigIRcoords = (zeros(size(ISCOORDS)));
tp@0: 	IRreflist = zeros(size(iv));
tp@0: 	
tp@0: 	% IMPROVE: Not necessary to calculate the IR coordinates for all
tp@0: 	% combinations since many are the same!
tp@0: 
tp@0:     % Now we make a temporary copy of POTENTIALISES which is displaced to
tp@0:     % the right since this will make it easier to align the postspec
tp@0:     % combos.
tp@0:     
tp@0:     PotentialISESshift = POTENTIALISES;
tp@0:     for ii = 2:specorder-1
tp@0:         iv = uint32(find(reflorder_select==ii));    
tp@0:         PotentialISESshift(masterivlistorig(iv),:) = [zeros(length(iv),specorder-ii) PotentialISESshift(masterivlistorig(iv),1:ii)];
tp@0:     end
tp@0:     
tp@0:     % Go through ii, which is the number of specular reflections after
tp@0:     % the diffraction.
tp@0: 
tp@0:     if ~isempty(nspecular_select)
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:             % We save the index vectors from the different orders so they can
tp@0:             % be reused in the next stage.
tp@0: 			
tp@0:              iv = uint32(find(nspecular_select == ii));
tp@0:              listselect = masterivlistorig(iv);
tp@0:              nlist = length(listselect);
tp@0:              eval(['iv',JJ(ii,1:JJnumbofchars(ii)),' = iv;']) 
tp@0:                  
tp@0:             if ii == 1
tp@0:                 bigIRcoords(listselect,:) = IRcoordsstart( PotentialISESshift(listselect,specorder),:);
tp@0:         
tp@0:             elseif ii == 2
tp@0:                 bigIRcoords(listselect,:) = EDB1findis(IRcoordsstart(PotentialISESshift(listselect,specorder),:),...
tp@0:                                                                PotentialISESshift(listselect,specorder-1),planeeqs,nlist,onesvec3);
tp@0: 	
tp@0:              elseif ii >= 3
tp@0:                 ivcombo = double(PotentialISESshift(listselect,specorder));
tp@0:                 for jj = 1:ii-2
tp@0:                     ivcombo = ivcombo + double(PotentialISESshift(listselect,specorder-jj))*ndecimaldivider^(jj);              
tp@0:                 end
tp@0:                 IRreflist = PointertoIRcombs(ivcombo);
tp@0:                 ivsimple    = find(IRreflist~=0);
tp@0:                 ivnotsimple = find(IRreflist==0);
tp@0:                 newIRcoordssimple    = EDB1findis(bigIRcoords(IRreflist(ivsimple),:),...
tp@0:                                                   PotentialISESshift(listselect(ivsimple),specorder-(ii-1)),planeeqs,size(ivsimple,1),onesvec3);
tp@0:                 newIRcoordsnotsimple = EDB1findis(IRcoordsstart(PotentialISESshift(listselect(ivnotsimple),specorder),:),...
tp@0:                                                   PotentialISESshift(listselect(ivnotsimple),specorder-1),planeeqs,size(ivnotsimple,1),onesvec3);
tp@0:                 for jj = 2:ii-1
tp@0:                     newIRcoordsnotsimple = EDB1findis(newIRcoordsnotsimple,...
tp@0:                                                   PotentialISESshift(listselect(ivnotsimple),specorder-jj),planeeqs,size(ivnotsimple,1),onesvec3);
tp@0:                 end
tp@0:                 newIRcoordstoadd = zeros(length(listselect),3);
tp@0:                 newIRcoordstoadd(ivsimple,:) = newIRcoordssimple;
tp@0:                 newIRcoordstoadd(ivnotsimple,:) = newIRcoordsnotsimple;
tp@0:                 bigIRcoords(listselect,:) = newIRcoordstoadd;
tp@0: 	
tp@0:             end
tp@0: 		
tp@0: 		end
tp@0:         clear nspecular_select
tp@0:     else
tp@0:         nspecular_select_not_empty = 0;
tp@0:     end
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: 
tp@0: %   ##################################################
tp@0: %   #                                                #
tp@0: %   #  S - (spec) - edge - spec - spec - R cases     #
tp@0: %   #                                                #
tp@0: %   #   (Pre- and) Postspec cases                    #
tp@0: %   #                                                #
tp@0: %   ##################################################
tp@0: %
tp@0: % Possible edges for S-E-spec-spec-R are seen (at least partly) by the receiver.
tp@0: %
tp@0: % The visibility must be checked, and also obstructions, but only
tp@0: % for the paths between the edge and the receiver
tp@0: 
tp@0: % The vector masterivlist will always refer to the original data vector
tp@0: % i.e. PotentialIR, OriginsfromIR, reflorderIR etc
tp@0: %
tp@0: % First we pick out those indices where there was a single diffraction, 
tp@0: % which wasn't last in the sequence 
tp@0: 
tp@0: if specorder > 1 & nspecular_select_not_empty == 1
tp@0: 
tp@0:     % The ii-value refers to the number of specular reflections after the
tp@0:     % diffraction
tp@0: 
tp@0: 	for ii = 1:specorder-1
tp@0:         
tp@0:         if SHOWTEXT >= 3
tp@0:             disp(['      Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the edge diff'])    
tp@0:         end
tp@0: 	
tp@0:         % The selection of cases with ii specular reflection after one
tp@0:         % diffraction have already been done, in the IR-building process.
tp@0: 	
tp@0:         eval(['masterivlist = masterivlistorig(iv',JJ(ii,1:JJnumbofchars(ii)),');'])
tp@0:         
tp@0:         if nedges < 256
tp@0:             possibleedges = uint8(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0:         elseif nedges < 65536
tp@0:             possibleedges = uint16(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0:         else
tp@0:             possibleedges = uint32(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0:         end
tp@0:         
tp@0:         % NB! possiblecombs will contain the specular sequences after the
tp@0:         % diffraction but they appear in reversed order!
tp@0:         possiblecombs = PotentialISESshift(masterivlist,specorder-ii+1:specorder);   
tp@0: 
tp@0:         % Specular combinations before the diffraction
tp@0:         nmaxprespecs = specorder-1-ii;
tp@0:         if nmaxprespecs == 0
tp@0:             if nedges+nplanes < 256
tp@0:                 prespeccombs = uint8(zeros(length(masterivlist),1));
tp@0:             elseif nedges+nplanes < 65536
tp@0:                 prespeccombs = uint16(zeros(length(masterivlist),1));
tp@0:             else
tp@0:                 prespeccombs = uint32(zeros(length(masterivlist),1));                
tp@0:             end
tp@0:             nmaxprespecs = 1;
tp@0:         else
tp@0:             prespeccombs = PotentialISESshift(masterivlist,1:nmaxprespecs);
tp@0:         end
tp@0:        % Visibility of the edge in each sequence, seen from the source
tp@0:         edgeweightsfromsou = ISESVISIBILITY(masterivlist);
tp@0:         
tp@0:         %--------------------------------------------------------------
tp@0:         % Expand to take the edge subdivisions into account - 
tp@0:         % treat all the edge subdivisions as separate edges for the
tp@0:         % visibility and obstruction test.
tp@0:         
tp@0:         nposs = length(masterivlist);
tp@0:         nposs = nposs*nedgesubs;
tp@0:        
tp@0:         expandedposscombs = repmat(possiblecombs,[1,nedgesubs]);
tp@0:         expandedposscombs = reshape(expandedposscombs.',ii,nposs).';
tp@0:   
tp@0:         expandedprespecs = repmat(prespeccombs,[1,nedgesubs]);
tp@0:         expandedprespecs = reshape(expandedprespecs.',nmaxprespecs,nposs).';
tp@0:                         
tp@0:         expandededgeweightsfromsou = repmat(edgeweightsfromsou,[1,nedgesubs]);
tp@0:         expandededgeweightsfromsou = reshape(expandededgeweightsfromsou.',1,nposs).';
tp@0:         
tp@0:         expandedpossedges = repmat(possibleedges,[1,nedgesubs]);
tp@0:         expandedpossedges = reshape(expandedpossedges.',1,nposs).';
tp@0: 
tp@0:         expandedmasterivlist = repmat(masterivlist,[1,nedgesubs]);
tp@0:         expandedmasterivlist = reshape(expandedmasterivlist.',1,nposs).';
tp@0: 
tp@0:         if SHOWTEXT >= 3
tp@0: 			disp(['         ',int2str(nposs),' Edge+IR segments found initially '])
tp@0:         end
tp@0: 	        
tp@0:         % Go through, iteratively, and check if the path from edge to R passes
tp@0:         % through all reflection planes along the way
tp@0:     
tp@0:        for jj = ii:-1:1
tp@0: 	
tp@0:             if nposs > 0
tp@0: 	
tp@0:                 if jj == ii
tp@0:                     fromcoords = full(bigIRcoords(expandedmasterivlist,:));
tp@0:                     [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0:                     tocoords = toedgecoords;
tp@0:                 else
tp@0:                     eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])    
tp@0:                     startlistvec = [1:length(expandedmasterivlist)];                    
tp@0:                     ivref = IRoriginsfrom(expandedmasterivlist);
tp@0:                     for kk = jj:ii-2
tp@0:                         ivnoIRexist = find(ivref==0);
tp@0:                         if ~isempty(ivnoIRexist)
tp@0:                             ivIRexist = startlistvec;
tp@0:                             ivIRexist(ivnoIRexist) = [];
tp@0:                             ivref(ivIRexist) = IRoriginsfrom(ivref(ivIRexist));
tp@0:                         else
tp@0:                             ivref = IRoriginsfrom(ivref);                            
tp@0:                         end
tp@0:                     end
tp@0:                     ivnoIRexist = find(ivref==0);                   
tp@0:                     if isempty(ivnoIRexist)
tp@0:                         fromcoords = full(bigIRcoords(ivref,:)); 
tp@0:                     else
tp@0:                         ivIRexist = startlistvec;
tp@0:                         ivIRexist(ivnoIRexist) = [];
tp@0:                         fromcoords = zeros(length(ivref),3);
tp@0:                         fromcoords(ivIRexist,:) = full(bigIRcoords(ivref(ivIRexist),:));
tp@0: 
tp@0:                         ISEScutout = PotentialISESshift(expandedmasterivlist(ivnoIRexist),specorder-ii+1:specorder);
tp@0:                         newIRcoords = EDB1findis(R,ISEScutout(:,ii),planeeqs,1,onesvec3);
tp@0:                         for kk = 2:jj
tp@0:                             newIRcoords = EDB1findis(newIRcoords,ISEScutout(:,ii-kk+1),planeeqs,size(newIRcoords,1),onesvec3);
tp@0:                         end
tp@0:                         fromcoords(ivnoIRexist,:) = newIRcoords;
tp@0:                         
tp@0:                     end
tp@0:                     
tp@0:                  end
tp@0:                 
tp@0:                 colno = ii-jj+1;
tp@0: 	
tp@0:                 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,colno),4),planenvecs(expandedposscombs(:,colno),:),minvals(expandedposscombs(:,colno),:),...
tp@0: 				    maxvals(expandedposscombs(:,colno),:),planecorners(expandedposscombs(:,colno),:),corners,ncornersperplanevec(expandedposscombs(:,colno)));
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:                 expandedmasterivlist  = expandedmasterivlist(hitplanes);
tp@0:                 expandedposscombs     = expandedposscombs(hitplanes,:);
tp@0:                 expandedpossedges     = expandedpossedges(hitplanes);
tp@0:                 expandedprespecs      = expandedprespecs(hitplanes,:);
tp@0:                 edgeweightlist        = edgeweightlist(hitplanes);
tp@0:                 expandededgeweightsfromsou = expandededgeweightsfromsou(hitplanes); 
tp@0:                 toedgecoords = toedgecoords(hitplanes,:);
tp@0: 	
tp@0:                 if jj < ii
tp@0:                     for kk = jj+1:ii
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(expandedmasterivlist);
tp@0:                 
tp@0:             end
tp@0:             if SHOWTEXT >= 3
tp@0: 			    disp(['         ',int2str(nposs),' Edge+IR segments survived the visibility test in refl plane ',JJ(jj,1:JJnumbofchars(jj))])               
tp@0:                 if SHOWTEXT >= 5
tp@0:                    POTENTIALISES(expandedmasterivlist,:)
tp@0:                 end
tp@0:             end
tp@0:             
tp@0:         end
tp@0:         
tp@0:         if obstructtestneeded & nposs > 0
tp@0:         
tp@0:             % Check obstructions for all the paths: edge -> plane1 -> plane2 -> ...
tp@0:             % -> edge   (S to edge is not needed!)
tp@0:             for jj = 1:ii+1
tp@0:                 if nposs > 0
tp@0:                     if jj==1
tp@0:                         fromcoords = R;    
tp@0:                         startplanes = [];    
tp@0:                     else
tp@0:                         startplanes = expandedposscombs(:,ii-jj+2);
tp@0:                         eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0:                     end
tp@0:                     if jj == ii+1
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(:,ii-jj+1);    
tp@0:                     end
tp@0:                                         
tp@0:                     [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0:                         planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0:                     if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0:                         disp('WARNING! An edgehit or cornerhit occurred during the obstruction test but this is not')
tp@0:                         disp('         handled correctly yet.')
tp@0:                     end
tp@0: 	
tp@0:                     if nobstructions > 0
tp@0:                         expandedmasterivlist  = expandedmasterivlist(nonobstructedpaths);
tp@0:                         expandedposscombs     = expandedposscombs(nonobstructedpaths,:);
tp@0:                         expandedpossedges     = expandedpossedges(nonobstructedpaths);
tp@0:                         expandedprespecs      = expandedprespecs(nonobstructedpaths,:);
tp@0:                         edgeweightlist        = edgeweightlist(nonobstructedpaths);
tp@0:                         expandededgeweightsfromsou = expandededgeweightsfromsou(nonobstructedpaths); 
tp@0:                         toedgecoords          = toedgecoords(nonobstructedpaths,:);
tp@0:                         nposs = length(expandedmasterivlist);
tp@0:                         
tp@0:                         for kk = 1:ii
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: 	
tp@0:                 end
tp@0:                 if SHOWTEXT >= 3
tp@0: 			       disp(['         ',int2str(nposs),' Edge+IR segments survived the obstruction test for path segment ',int2str(jj)])
tp@0:                     if SHOWTEXT >= 5
tp@0:                         POTENTIALISES(expandedmasterivlist,:)
tp@0:                     end
tp@0:                 end
tp@0:            end      
tp@0:  
tp@0:         end
tp@0:         
tp@0:         if nposs > 0
tp@0: 
tp@0:             % Visibility of edge segments is given by bitand-ing the
tp@0:             % edge visibility from the (image) source and from the receiver
tp@0: 
tp@0:             edgeweightlist = bitand(edgeweightlist,expandededgeweightsfromsou);
tp@0:             iv = find(edgeweightlist==0);
tp@0: 
tp@0:             if ~isempty(iv)
tp@0:                 expandedpossedges(iv) = [];
tp@0:                 expandedposscombs(iv,:) = [];
tp@0:                 expandedprespecs(iv,:) = [];
tp@0:                 edgeweightlist(iv) = [];
tp@0:                 expandedmasterivlist(iv) = [];
tp@0:                 nposs = length(edgeweightlist);
tp@0:                 
tp@0:             end
tp@0:              
tp@0:             % Add the newly found combinations to the outdata list
tp@0:             
tp@0:             [nnew,slask] = size(expandedprespecs);
tp@0:             if nnew == 1
tp@0:                 colstoclear = find( expandedprespecs==0 );
tp@0:             else
tp@0:                 colstoclear = find(sum(expandedprespecs>0)==0);
tp@0:             end
tp@0:             if ~isempty(colstoclear)
tp@0:                 expandedprespecs(:,colstoclear) = [];    
tp@0:             end
tp@0: 
tp@0:             [slask,ncolsexist] = size(prespeclist);
tp@0:             [slask,ncolsnew] = size(expandedprespecs);
tp@0:             ncolstoadd = ncolsexist - ncolsnew;
tp@0:             if ncolstoadd > 0
tp@0:                 prespeclist =  [prespeclist; [expandedprespecs zeros(nposs,ncolstoadd)]];
tp@0:             else
tp@0:                 prespeclist =  [prespeclist; expandedprespecs];                
tp@0:             end
tp@0:             edgedifflist = [edgedifflist;expandedpossedges];
tp@0:             postspeclist = [postspeclist;[expandedposscombs zeros(nposs,specorder-1-ii)]];
tp@0:             bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0:             validIRcoords = [validIRcoords;bigIRcoords(expandedmasterivlist,:)];
tp@0:             % NB! It is correct below that the indices for the ISCOORDS should be
tp@0:             % ORIGINSFROM(masterivlist), rather than masterivlist.
tp@0:             % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0:             % spec-spec-...-diff combinations and then
tp@0:             % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0:             % ends with a diff has no image source. 
tp@0:             % If we have a spec-spec-...-diff-spec comb., then we must
tp@0:             % use ORIGINSFROM(ORIGINSFROM(masterivlist)) etc.
tp@0:             ivref = ORIGINSFROM(expandedmasterivlist);
tp@0:             for kk = 1:ii
tp@0:                 ivref = ORIGINSFROM(ivref);
tp@0:             end
tp@0:             validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0: 
tp@0:         end
tp@0: 	end
tp@0: 
tp@0: %    clear PotentialISESshift IScoords bigIRcoords
tp@0:     clear PotentialISESshift bigIRcoords
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #
tp@0: %   #   Pack the edge segments together because every little edge segment
tp@0: %   #   is present as a separate edge
tp@0: %   #   This can be done for all combinations at once.
tp@0: %   #
tp@0: %   #######################################################################
tp@0: 
tp@0: test = [prespeclist edgedifflist postspeclist];
tp@0: 
tp@0: if ~isempty(test)
tp@0: 
tp@0: 	ncombs = length(edgedifflist);
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:         bigedgeweightlist(ivremove+1) = double(bigedgeweightlist(ivremove+1)) + double(bigedgeweightlist(ivremove));
tp@0:         bigedgeweightlist(ivremove) = [];
tp@0:         edgedifflist(ivremove) = [];
tp@0:         prespeclist(ivremove,:) = [];
tp@0:         postspeclist(ivremove,:) = [];
tp@0:         validIScoords(ivremove,:) = [];
tp@0:         validIRcoords(ivremove,:) = [];
tp@0: 	
tp@0:         test = [prespeclist edgedifflist postspeclist];
tp@0:         ncombs = length(edgedifflist);
tp@0:         dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
tp@0:         ivremove = find(dtest==1);
tp@0: 	end
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #
tp@0: %   #   The weights of the visible edge segments should be
tp@0: %   #   translated into start and end points, together with the visibility
tp@0: %   #   weights from the receiver.
tp@0: %   #   This can be done for all combinations at once.
tp@0: %   #
tp@0: %   #######################################################################
tp@0: 
tp@0: % As a start we set the start and end values to 0 and 1, i.e. assuming full
tp@0: % visibility.
tp@0: 
tp@0: if specorder >= 1 & ~isempty(edgedifflist)
tp@0: 	startandendpoints = [startandendpoints;[zeros(length(edgedifflist),1) ones(length(edgedifflist),1)]];
tp@0: 	
tp@0:     totvisibility = bigedgeweightlist;
tp@0: 
tp@0: 	ivnovisibility = find(totvisibility==0);
tp@0: 	if ~isempty(ivnovisibility)
tp@0:         disp('Clearing some cases.')
tp@0:         edgedifflist(ivnovisibility) = [];
tp@0:         prespeclist(ivnovisibility,:) = [];
tp@0:         postspeclist(ivnovisibility,:) = [];
tp@0:         bigedgeweightlist(ivnovisibility) = [];
tp@0:         totvisibility(ivnovisibility) = [];
tp@0:         validIScoords(ivnovisibility,:) = [];
tp@0:         validIRcoords(ivnovisibility,:) = [];
tp@0: 	end    
tp@0: 	
tp@0: 	ivtoaccountfor = [1:length(totvisibility)].';
tp@0: 	
tp@0: 	% If there are edges with full visibility we don't need to do anything to
tp@0: 	% the start and end values, but we remove the edges from the list to
tp@0: 	% process.
tp@0: 	
tp@0: 	ivwholeedges = find( totvisibility == maxvisibilityvalue);
tp@0: 	if ~isempty(ivwholeedges)
tp@0:         ivtoaccountfor(ivwholeedges) = [];
tp@0: 	end
tp@0: 	
tp@0: 	if ~isempty(ivtoaccountfor)
tp@0:         ncombs = length(ivtoaccountfor);
tp@0:         bitpattern = zeros(ncombs,nedgesubs);
tp@0:         for ii=1:nedgesubs
tp@0:             bitpattern(:,ii) = bitget(totvisibility(ivtoaccountfor),ii); 
tp@0:         end
tp@0:         dbit1 = diff([zeros(ncombs,1) bitpattern].').';
tp@0:         dbit2 = [dbit1 -bitpattern(:,nedgesubs)]; 
tp@0:         
tp@0:         nsegments = ceil((sum(abs(dbit1.')).')/2);
tp@0: 	
tp@0:         ivonesegments = find(nsegments==1);
tp@0:         if ~isempty(ivonesegments)
tp@0: 
tp@0:             nonesegments = length(ivonesegments);
tp@0:             multvec = 2.^[0:nedgesubs];
tp@0:             segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0:             segendpos   = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0: 	
tp@0:             ivmodify = find(segstartpos==1);
tp@0:             segstartpos(ivmodify) = ones(size(ivmodify))*1.5;
tp@0:             ivmodify = find(segendpos>nedgesubs);
tp@0:             segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5);
tp@0: 	
tp@0:             startandendpoints(ivtoaccountfor(ivonesegments),1) = (segstartpos-1.5)/(nedgesubs-1);
tp@0:             startandendpoints(ivtoaccountfor(ivonesegments),2) = (segendpos-1.5)/(nedgesubs-1);
tp@0:                     
tp@0:         end    
tp@0:         
tp@0:         % If we have some two-or-more-subsegments cases, they will be
tp@0:         % discovered by the if-condition below
tp@0:         
tp@0:         if length(ivonesegments) < ncombs
tp@0:             for nsegmentstocheck = 2:ceil(nedgesubs/2)
tp@0:      
tp@0:                 ivNsegments = find(nsegments==nsegmentstocheck);
tp@0:                 if ~isempty(ivNsegments)
tp@0:                     [n1,n2] = size(startandendpoints);
tp@0:                     if n2 < 2*nsegmentstocheck
tp@0:                         startandendpoints = [startandendpoints zeros(n1,2*nsegmentstocheck-n2)];    
tp@0:                     end
tp@0:                     for jj = 1:length(ivNsegments)
tp@0:                         ivstartbits = find(dbit2(ivNsegments(jj),:) == 1);
tp@0:                         ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits;
tp@0:                         ivendbits = find(dbit2(ivNsegments(jj),:) == -1);
tp@0:                         ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits;
tp@0:                         
tp@0:                         for kk = 1:nsegmentstocheck,                                                        
tp@0:                             startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*2+1) = (ivstartbits(kk)-1.5)/(nedgesubs-1);
tp@0:                             startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*2+2) = (ivendbits(kk)-1.5)/(nedgesubs-1);
tp@0:                         end
tp@0:                     end                
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:         end        
tp@0: 	end
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #
tp@0: %   #   Construct a list guide, which will tell which rows have only
tp@0: %   #   d, which rows have sd etc
tp@0: %   #   
tp@0: %   #######################################################################
tp@0: 
tp@0: [n1,n2] = size(prespeclist);
tp@0: if n1 ~= 0
tp@0: 	if n2 == 0
tp@0:         prespeclist = zeros(n1,1);    
tp@0: 	end
tp@0: 	
tp@0: 	if n2 > 1
tp@0:         nprespecs  = sum(prespeclist.' > 0).';
tp@0: 	else
tp@0:         nprespecs = (prespeclist>0);
tp@0: 	end
tp@0: 	[n1,n2] = size(postspeclist);
tp@0: 	if n2 > 1
tp@0:         npostspecs  = sum(postspeclist.' > 0).';
tp@0: 	else
tp@0:         npostspecs = (postspeclist>0);    
tp@0: 	end
tp@0: 	
tp@0: 	listofdiffcol   = 1 + nprespecs;
tp@0: 	listofreflorder = listofdiffcol + npostspecs;
tp@0: 	
tp@0: 	[B,ivec,jvec] = unique([listofreflorder listofdiffcol],'rows');
tp@0: 	nuniquecombs = length(ivec);
tp@0: 	ntotcombs = length(jvec);
tp@0: 	
tp@0: 	listguide = zeros(nuniquecombs,3);
tp@0: 	listoforders = zeros(nuniquecombs,2);
tp@0: 	sortvec = zeros(ntotcombs,1);
tp@0: 	for ii = 1:length(ivec)
tp@0:         ivfindcombs = find(jvec==ii);
tp@0:         listguide(ii,1) = length(ivfindcombs);
tp@0:         if ii > 1
tp@0:             listguide(ii,2) = listguide(ii-1,3)+1;
tp@0:         else
tp@0:             listguide(ii,2) = 1;    
tp@0:         end
tp@0:         listguide(ii,3) = listguide(ii,2)+listguide(ii,1)-1;
tp@0:         listoforders(ii,:) = [B(ii,1) B(ii,2)]; 
tp@0: 	
tp@0:         sortvec(listguide(ii,2):listguide(ii,3)) = ivfindcombs;
tp@0:         
tp@0: 	end
tp@0: 	
tp@0: 	prespeclist = prespeclist(sortvec,:);
tp@0: 	postspeclist = postspeclist(sortvec,:);
tp@0: 	validIScoords = validIScoords(sortvec,:);
tp@0: 	validIRcoords = validIRcoords(sortvec,:);
tp@0: 	edgedifflist = edgedifflist(sortvec,:);
tp@0: 	startandendpoints = startandendpoints(sortvec,:);
tp@0: 	bigedgeweightlist = bigedgeweightlist(sortvec,:);
tp@0:     
tp@0:     % The prespeclist needs to be shifted back to the left because it is
tp@0:     % "right-aligned"
tp@0:     
tp@0:     [ncombs,ncols] = size(prespeclist);
tp@0:     if ncols > 1
tp@0:         if ncombs > 1
tp@0:             iv = find( ( prespeclist(:,1)==0 ) & ( prespeclist(:,2)~=0 ) );
tp@0:             if ~isempty(iv)
tp@0:                prespeclist(iv,1:ncols-1) = prespeclist(iv,2:ncols);
tp@0:                prespeclist(iv,ncols) = 0;
tp@0:             end
tp@0: 
tp@0:             for jj = 3:ncols
tp@0:                 iv = find( ( sum(prespeclist(:,1:jj-1).').' == 0 ) & ( prespeclist(:,jj)~=0 ) );
tp@0:                 if ~isempty(iv)
tp@0:                    prespeclist(iv,1:ncols-(jj-1)) = prespeclist(iv,jj:ncols);
tp@0:                    prespeclist(iv,ncols-(jj-2):ncols) = 0;
tp@0:                 end
tp@0:             end    
tp@0:         else
tp@0:            if prespeclist(1,1) == 0,               
tp@0:                 ninitialzeros =  find(diff(cumsum(prespeclist)));
tp@0:                 if ~isempty(ninitialzeros)
tp@0:                     ninitialzeros = ninitialzeros(1);
tp@0:                     prespeclist(1:ncols-ninitialzeros) = prespeclist(ninitialzeros+1:ncols);
tp@0:                     prespeclist(ncols-ninitialzeros+1:ncols) = 0;
tp@0:                 end
tp@0:            end
tp@0:         end
tp@0:     end
tp@0:     
tp@0: else
tp@0:     prespeclist = [];
tp@0:     postspeclist = [];
tp@0:     validIScoords = [];
tp@0:     validIRcoords = [];
tp@0:     edgedifflist = [];
tp@0:     startandendpoints = [];
tp@0:     bigedgeweightlist = [];
tp@0:     listguide = [];
tp@0:     
tp@0: end