tp@0: function [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,listguide,...
tp@0:     listoforders] = EDB1diffNISES(eddatafile,S,R,...
tp@0:     ivNdiffmatrix,lengthNdiffmatrix,Ndifforder,specorder,visplanesfromR,vispartedgesfromS,...
tp@0:     vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlistin,validEDIRcoords)
tp@0: % EDB1diffNISES - Gives a list of paths that includes an N-order diffraction path.
tp@0: % Gives the list of visible N-diffraction paths, possibly with specular reflections before and after.
tp@0: %
tp@0: % Input parameters:
tp@0: %    eddatafile,S,R,ivsinglediff,...
tp@0: %    singlediffcol,startindicessinglediff,endindicessinglediff,specorder,visplanesfromR,...
tp@0: %    vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,PointertoIRcombs,IRoriginsfrom;
tp@0: %
tp@0: % GLOBAL parameters:
tp@0: %   SHOWTEXT JJ JJnumbofchars   See EDB1mainISES
tp@0: %   POTENTIALISES,ISCOORDS,ORIGINSFROM,ISESVISIBILITY,REFLORDER   See EDB1findISEStree
tp@0: %
tp@0: % Output parameters:
tp@0: %   edgedifflist        List [ncombs,N] of the edge numbers involved in each
tp@0: %                       spec-diff-...-diff-spec combination.
tp@0: %   startandendpoints   Matrix [ncombs,4] of the relative start and end
tp@0: %                       points of each edge. The values, [0,1], indicate
tp@0: %                       which part of the two edges that are visible.
tp@0: %   prespeclist         Matrix [ncombs,specorder-N] of the specular
tp@0: %                       reflections that precede every diffraction.
tp@0: %   postspeclist        Matrix [ncombs,specorder-N] 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-...-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-...-diff-spec comb.
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) 20031018
tp@0: %
tp@0: % [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,...
tp@0: %   listguide] = EDB1diff3ISES(eddatafile,S,R,...
tp@0: %   ivNdiffmatrix,lengthNdiffmatrix,Ndifforder,specorder,visplanesfromR,vispartedgesfromS,...
tp@0: %   vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlist,validEDIRcoords);
tp@0: 
tp@0: global SHOWTEXT JJ JJnumbofchars
tp@0: global POTENTIALISES ISCOORDS ORIGINSFROM ISESVISIBILITY REFLORDER
tp@0: 
tp@0: eval(['load ',eddatafile])
tp@0: 
tp@0: [nedges,slask] = size(planesatedge);
tp@0: [nplanes,slask] = size(planecorners);
tp@0: multfac = 10^(ceil(log10(max([nedges nplanes]))));
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,specorder-Ndifforder);
tp@0: zerosvec2 = zeros(1,3);
tp@0: listguide = zeros(specorder*2-1,3);
tp@0: 
tp@0: obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0: onesvec = ones(1,nedgesubs);
tp@0: onesvec3 = ones(1,3);
tp@0: 
tp@0: npostspecs = sum(double(postspeclistin.'>0)).';
tp@0: 
tp@0: %   ###########################################
tp@0: %   #                                         #
tp@0: %   #   S - edge - edge - edge -R cases       #
tp@0: %   #                                         #
tp@0: %   ###########################################
tp@0: %
tp@0: % Possible edges for S-E-E-E-R are seen (at least partly) by the source and by the receiver.
tp@0: %
tp@0: % The visibility by the source is taken care of by the findISEStree
tp@0: % so we must check which ones are visible by the receiver.
tp@0: % Also, the active edge segment(s) must be selected but this is done
tp@0: % further down together with the S-spec-spec-edge-R cases
tp@0: 
tp@0: ivNdiff = ivNdiffmatrix(1:lengthNdiffmatrix(Ndifforder),Ndifforder);
tp@0: ivsinglediff = ivNdiffmatrix(1:lengthNdiffmatrix(1),1);
tp@0: 
tp@0: ivSEEER = ivNdiff(find(REFLORDER(ivNdiff)==Ndifforder));
tp@0: 
tp@0: possibleedgecombs = double(POTENTIALISES(ivSEEER,1:Ndifforder))-nplanes;
tp@0: ivnotvisiblefromr = find(vispartedgesfromR(possibleedgecombs(:,Ndifforder))==0);
tp@0: if ~isempty(ivnotvisiblefromr)
tp@0:     possibleedgecombs(ivnotvisiblefromr,:) = [];
tp@0: end
tp@0: 
tp@0: edgedifflist      = [edgedifflist;possibleedgecombs];
tp@0: bigedgeweightlist = [bigedgeweightlist;[vispartedgesfromS(edgedifflist(:,1)) vispartedgesfromR(edgedifflist(:,Ndifforder))]];
tp@0: 
tp@0: [nedgesadded,slask] = size(edgedifflist);
tp@0: zerosvec3 = zeros(nedgesadded,1);
tp@0: ndiffonly = nedgesadded;
tp@0: 
tp@0: prespeclist =  [prespeclist;zerosvec3(:,ones(1,max(specorder-Ndifforder,1)))];
tp@0: postspeclist = [postspeclist;zerosvec3(:,ones(1,max(specorder-Ndifforder,1)))];
tp@0: validIScoords = [validIScoords;S(ones(nedgesadded,1),:)];
tp@0: validIRcoords = [validIRcoords;R(ones(nedgesadded,1),:)];
tp@0: 
tp@0: if SHOWTEXT >= 3
tp@0: 	disp(['         ',int2str(nedgesadded),' Nth diff valid'])
tp@0: end
tp@0: 
tp@0: %   ###########################################
tp@0: %   ###########################################
tp@0: %   ###########################################
tp@0: %   ###########################################
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 the source-to-edge paths
tp@0: % were checked in the ISEStree, and the visibility from the receiver also
tp@0: % has been checked.
tp@0: 
tp@0: % The vector ivmultidiff will always refer to the original data vector
tp@0: % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc
tp@0: %
tp@0: % We should remove the combinations which involve an edge that the
tp@0: % receiver can not see, but since the edge number is in different columns
tp@0: % we do that in the for loop.
tp@0: 
tp@0: % The ii-loop will go through: spec-diff, spec-spec-diff
tp@0: % spec-spec-spec-diff etc
tp@0: 
tp@0: for ii = 1:specorder-Ndifforder
tp@0:     
tp@0:     if SHOWTEXT >= 3
tp@0:         disp(['      Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before the  Nth edge diff'])    
tp@0:     end
tp@0: 
tp@0:     % Select the combinations where the reflection order == ii+Ndifforder
tp@0:     % (which means ii specular reflections in addition to the diffraction)
tp@0:     % and where the last Ndifforder columns contain edges.
tp@0:     
tp@0:     iv = find((REFLORDER(ivNdiff)==ii+Ndifforder) & prod(   double( POTENTIALISES(ivNdiff,ii+1:ii+Ndifforder)>nplanes).'  ).'   );
tp@0:     masterivlist = ivNdiff(iv);
tp@0:     possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0:     
tp@0:     % Keep only combinations for which the receiver can see the edge
tp@0:     
tp@0:     ivnotvisiblefromr = find(vispartedgesfromR(possibleedgecombs(:,Ndifforder))==0);
tp@0:     if ~isempty(ivnotvisiblefromr)
tp@0:         masterivlist(ivnotvisiblefromr) = [];
tp@0:         possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0:     end        
tp@0: 
tp@0:     possibleprespecs = POTENTIALISES(masterivlist,1:ii);
tp@0:     
tp@0:     reftoIScoords = ORIGINSFROM(masterivlist);
tp@0:     edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgecombs(:,Ndifforder))];
tp@0:         
tp@0:     nposs = length(masterivlist);
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0:  		disp(['         ',int2str(nposs),' IS+edge multiples valid'])
tp@0:     end
tp@0: 
tp@0:     edgedifflist = [edgedifflist;possibleedgecombs];
tp@0:     prespeclist = [prespeclist;[possibleprespecs zeros(nposs,specorder-Ndifforder-ii)]];        
tp@0:     postspeclist = [postspeclist;zeros(nposs,specorder-Ndifforder)];
tp@0:     bigedgeweightlist = [bigedgeweightlist;edgeweightlist];    
tp@0:     % The ref. to ISCOORDS must be nested the correct number of time to
tp@0:     % come back to the image source coords of the last spec. refl.
tp@0:     ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist)));
tp@0:     for jj = 4:Ndifforder
tp@0:         ivref = ORIGINSFROM(ivref);    
tp@0:     end
tp@0:     validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0:     validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0:     
tp@0: end
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: 
tp@0: %   #############################################
tp@0: %   #                                           #
tp@0: %   #         S - edge - edge -spec - R cases   #
tp@0: %   #                                           #
tp@0: %   #         Postspec cases                    #
tp@0: %   #                                           #
tp@0: %   #############################################
tp@0: %
tp@0: % Possible edges for S-E-E-spec-spec-R are seen (at least partly) by the receiver.
tp@0: %
tp@0: % For some of the combos, we have the IR coords and the edge visibility
tp@0: % from the single-diffraction run. For potential combos that were not
tp@0: % included there, they are either invisible/obstructed or were not tested.
tp@0: % We can figure out which ones were tested because they can be found in the
tp@0: % POTENTIALISES under edge-spec-spec but not in the final valid list.
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, 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 first in the sequence
tp@0: % of reflections.
tp@0: 
tp@0: for ii = 1:specorder-Ndifforder
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0:         disp(['      Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the multiple edge diff'])    
tp@0:     end
tp@0: 
tp@0:     % Select the combinations where the reflection order == ii+Ndifforder
tp@0:     % (which means ii specular reflections in addition to the diffraction)
tp@0:     % and where the first Ndifforder columns contain edges.
tp@0:    
tp@0:     iv = find((REFLORDER(ivNdiff)==ii+Ndifforder) & prod(   double( POTENTIALISES(ivNdiff,1:Ndifforder)>nplanes).'  ).'   );
tp@0:     masterivlist = ivNdiff(iv);
tp@0:     possibleedgecombs = double(POTENTIALISES(masterivlist,1:Ndifforder)) - nplanes;
tp@0:     possiblepostspecs = POTENTIALISES(masterivlist,Ndifforder+1:specorder);
tp@0:     possibleweights = ISESVISIBILITY(masterivlist);
tp@0:     
tp@0:     % Compare with those that have already been found OK
tp@0:     ivOK = find(npostspecs==ii);
tp@0:     if ~isempty(ivOK)
tp@0:         patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:ii)];
tp@0:     else
tp@0:         patternOK = [];    
tp@0:     end
tp@0:     % Find out which ones have been checked and found invisible/obstructed
tp@0:     ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == ii+1));
tp@0:     patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+ii))];
tp@0:     if ~isempty(patternOK)
tp@0:         patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0:     else
tp@0:         patternNOTOK = patternALL; 
tp@0:     end
tp@0:     patterntocompare = [possibleedgecombs(:,Ndifforder) possiblepostspecs(:,1:ii)];
tp@0:     
tp@0:     ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0:     masterivlist(ivtocancel) = [];
tp@0:     possibleedgecombs(ivtocancel,:) = [];
tp@0:     possiblepostspecs(ivtocancel,:) = [];
tp@0:     possibleweights(ivtocancel,:) = [];
tp@0:     patterntocompare(ivtocancel,:) = [];
tp@0: 
tp@0:     [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0:     ivmustbechecked = find(ivcompletelyOK==0);
tp@0:     ivcompletelyOK = find(ivcompletelyOK);
tp@0:     if ~isempty(ivmustbechecked)
tp@0:         masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0:         ntocheckmore = length(masterlisttocheckmore);
tp@0:         
tp@0:         %----------------------------------------------
tp@0:         % Must carry out a visibility and obstruction check for the special
tp@0:         % combinations here.
tp@0:         % These combinations have a postspec-combination that hasn't been
tp@0:         % encountered in the single diffraction cases, so no visibility
tp@0:         % test has been made for these.
tp@0:         
tp@0:         lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,Ndifforder))-nplanes;
tp@0:         newIRcoords = R;
tp@0:         reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0:         for jj = 1:ii
tp@0:             reflplanes = POTENTIALISES(masterlisttocheckmore,Ndifforder+1+ii-jj);
tp@0:             reflplanesexpand(:,jj) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0:             newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,size(newIRcoords,1),onesvec3);
tp@0:             newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0:             eval(['newIRcoords',JJ(jj,1:JJnumbofchars(jj)),' = newIRcoordsexpand;'])                    
tp@0:         end
tp@0:         [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0:         tocoords = toedgecoords;
tp@0:         lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0:         lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0:         masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0:         masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0:         
tp@0:         ntocheckmore = length(masterlisttocheckmore);
tp@0:         if SHOWTEXT >= 3
tp@0:             disp(['         ',int2str(ntocheckmore),' special N-edge - IR combinations to check'])    
tp@0:         end
tp@0:         
tp@0:         for jj = ii:-1:1
tp@0:             if length(masterlisttocheckmore) > 0
tp@0:                 eval(['fromcoords = newIRcoords',JJ(jj,1:JJnumbofchars(jj)),';']);
tp@0:                 if jj < ii
tp@0:                     eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])    
tp@0:                 end
tp@0:                 
tp@0:                 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,jj),4),planenvecs(reflplanesexpand(:,jj),:),minvals(reflplanesexpand(:,jj),:),...
tp@0: 				    maxvals(reflplanesexpand(:,jj),:),planecorners(reflplanesexpand(:,jj),:),corners,ncornersperplanevec(reflplanesexpand(:,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('         implemented.')
tp@0:                 end
tp@0:                 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0: 
tp@0:                 masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0:                 edgeweightlist = edgeweightlist(hitplanes);
tp@0:                 lastedgenumbers = lastedgenumbers(hitplanes);
tp@0:                 reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0:                 toedgecoords = toedgecoords(hitplanes,:);
tp@0: 
tp@0:                 for kk = 1:ii
tp@0:                     eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0:                     if kk > jj
tp@0:                         eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0:                     end
tp@0:                 end
tp@0:                 ntocheckmore = length(masterlisttocheckmore);
tp@0:             
tp@0:             end
tp@0:             
tp@0:             if SHOWTEXT >= 3
tp@0:                 disp(['         ',int2str(ntocheckmore),' of the special N-edge - IR combinations survived the visibility test in refl plane ',int2str(jj)])
tp@0:             end
tp@0:         end
tp@0: 
tp@0:         if obstructtestneeded & ntocheckmore > 0
tp@0: 
tp@0:             for jj = 1:ii+1
tp@0:                 if ntocheckmore > 0
tp@0:                     if jj == 1
tp@0:                         fromcoords = R;
tp@0:                         startplanes = [];
tp@0:                     else
tp@0:                         eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';']) 
tp@0:                         startplanes = reflplanesexpand(:,jj-1);
tp@0:                     end
tp@0:                     if jj == ii+1,                    
tp@0:                         tocoords = toedgecoords;
tp@0:                         endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0:                     else
tp@0:                         eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])    
tp@0:                         endplanes = reflplanesexpand(:,jj);    
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('         implemented.')
tp@0:                     end
tp@0: 
tp@0:                     if nobstructions > 0
tp@0:                         masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0:                         edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0:                         lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0:                         reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0:                         toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0:                         for kk = 1:ii
tp@0:                             eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0:                             eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0:                         end
tp@0:                         
tp@0:                     end
tp@0:                     ntocheckmore = length(masterlisttocheckmore);
tp@0:                     
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:             if SHOWTEXT >= 3
tp@0:                 disp(['         ',int2str(ntocheckmore),' of the special N-edge - IR combinations survived the obstruction test'])
tp@0:             end
tp@0:             
tp@0:         end
tp@0: 
tp@0:         % Add the found special combinations to the outdata list
tp@0:     
tp@0:         edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,1:Ndifforder))-nplanes];
tp@0:         prespeclist = [prespeclist;zerosvec1(ones(ntocheckmore,1),:)];
tp@0:         postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0:         bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];    
tp@0:         eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(ii,1:JJnumbofchars(ii)),'];']);
tp@0:         validIScoords = [validIScoords;S(ones(ntocheckmore,1),:)];
tp@0:         
tp@0:     end
tp@0:     
tp@0:     masterivlist = masterivlist(ivcompletelyOK);    
tp@0:     possibleedgecombs = possibleedgecombs(ivcompletelyOK,:);
tp@0:     possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0:     possibleweights = possibleweights(ivcompletelyOK,:);
tp@0:         
tp@0:     nposs = length(ivcompletelyOK);
tp@0: 
tp@0:     if SHOWTEXT >= 3
tp@0: 		disp(['         ',int2str(nposs),' Edge,Nth-order + IR segments (non-special) survived the obstruction test'])
tp@0:     end
tp@0:     
tp@0:     % Add the found "standard" combinations to the outdata list
tp@0: 
tp@0:     edgedifflist = [edgedifflist;possibleedgecombs];
tp@0:     prespeclist = [prespeclist;zeros(nposs,specorder-Ndifforder)];
tp@0:     postspeclist = [postspeclist;possiblepostspecs ];
tp@0:     bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];    
tp@0:     validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0:     validIScoords = [validIScoords;S(ones(nposs,1),:)];
tp@0:     
tp@0: end
tp@0: 
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: %   #######################################################################
tp@0: 
tp@0: %   #############################################
tp@0: %   #                                           #
tp@0: %   #         S - spec - edge - spec - R cases  #
tp@0: %   #                                           #s
tp@0: %   #         Pre and postspec cases            #
tp@0: %   #                                           #
tp@0: %   #############################################
tp@0: 
tp@0: % The ii- and jj-loops will go through all combinations of spec before and after
tp@0: % the Ndifforder diffractions.
tp@0: %
tp@0: % Unlike before, we will look in the list of already found combinations
tp@0: % (edgedifflist etc)
tp@0: 
tp@0: 
tp@0: for ii = 1:specorder-Ndifforder-1
tp@0: 
tp@0:     for jj = 1:specorder-ii-Ndifforder
tp@0:     
tp@0:         if SHOWTEXT >= 3
tp@0:             disp(['      Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before and ',JJ(jj,1:JJnumbofchars(jj)),' spec refl after the N-order edge diff'])    
tp@0:         end
tp@0: 	
tp@0:         % Select the combinations where the reflection order ==
tp@0:         % ii+jj+Ndifforder
tp@0:         % (which means ii+jj specular reflections in addition to the
tp@0:         % diffraction), and where columns ii+1:ii+Ndifforder of POTENTIALISES
tp@0:         % contain edges.
tp@0: 	
tp@0:         iv = find((REFLORDER(ivNdiff) == ii+jj+Ndifforder) & prod(   double( POTENTIALISES(ivNdiff,ii+1:ii+Ndifforder)>nplanes).'  ).'      );
tp@0:         masterivlist = ivNdiff(iv);
tp@0:         possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0:         possibleprespecs  = POTENTIALISES(masterivlist,1:ii);
tp@0:         possiblepostspecs = POTENTIALISES(masterivlist,ii+Ndifforder+1:ii+Ndifforder+jj);
tp@0:         possibleweights   = ISESVISIBILITY(masterivlist);
tp@0: 
tp@0:         % Compare with those that have already been found OK
tp@0:         ivOK = find(npostspecs==jj);
tp@0:         if ~isempty(ivOK)
tp@0:             patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:jj)];
tp@0:         else
tp@0:             patternOK = [];
tp@0:         end
tp@0:         
tp@0:         % Find out which ones have been checked and found invisible/obstructed
tp@0:         ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == jj+1));
tp@0:         patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+jj))];
tp@0:         if (   (~isempty(patternOK)) &  (~isempty(patternALL)) )
tp@0:             patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0:         else
tp@0:             if (   (~isempty(patternOK)) )
tp@0:                 patternNOTOK = patternOK;
tp@0:             else
tp@0:                 patternNOTOK = patternALL;
tp@0:             end
tp@0:         end
tp@0:         patterntocompare = [possibleedgecombs(:,Ndifforder) possiblepostspecs(:,1:jj)];
tp@0: 
tp@0:         ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0:         masterivlist(ivtocancel) = [];
tp@0:         possibleedgecombs(ivtocancel,:) = [];
tp@0:         possibleprespecs(ivtocancel,:) = [];
tp@0:         possiblepostspecs(ivtocancel,:) = [];
tp@0:         possibleweights(ivtocancel,:) = [];
tp@0:         patterntocompare(ivtocancel,:) = [];
tp@0: 	
tp@0:         [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0:         ivmustbechecked = find(ivcompletelyOK==0);
tp@0:         ivcompletelyOK = find(ivcompletelyOK);
tp@0:         if ~isempty(ivmustbechecked)
tp@0:             masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0:             ntocheckmore = length(masterlisttocheckmore);
tp@0:         
tp@0:             %----------------------------------------------
tp@0:             % Must carry out a visibility and obstruction check for the special
tp@0:             % combinations here.
tp@0: 
tp@0:             lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,ii+Ndifforder))-nplanes;
tp@0:             newIRcoords = R;
tp@0:             reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0:             for kk = 1:jj
tp@0:                 reflplanes = POTENTIALISES(masterlisttocheckmore,Ndifforder+1+ii+jj-kk);
tp@0:                 reflplanesexpand(:,kk) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0:                 newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,size(newIRcoords,1),onesvec3);
tp@0:                 newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0:                 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoordsexpand;'])                    
tp@0:             end
tp@0:             [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0:             tocoords = toedgecoords;
tp@0:             lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0:             lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0:             masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0:             masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0:             
tp@0:             ntocheckmore = length(masterlisttocheckmore);
tp@0:             if SHOWTEXT >= 3
tp@0:                 disp(['         ',int2str(ntocheckmore),' special IS - N-edge - IR combinations to check'])    
tp@0:             end
tp@0:             
tp@0:             for kk = jj:-1:1
tp@0:                 if length(masterlisttocheckmore) > 0
tp@0:                     eval(['fromcoords = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),';']);
tp@0:                     if kk < jj
tp@0:                         eval(['tocoords = reflpoints',JJ(kk+1,1:JJnumbofchars(kk+1)),';'])    
tp@0:                     end
tp@0:                     
tp@0:                     [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,kk),4),planenvecs(reflplanesexpand(:,kk),:),minvals(reflplanesexpand(:,kk),:),...
tp@0: 					    maxvals(reflplanesexpand(:,kk),:),planecorners(reflplanesexpand(:,kk),:),corners,ncornersperplanevec(reflplanesexpand(:,kk)));
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('         implemented.')
tp@0:                     end
tp@0:                     eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints;'])
tp@0: 
tp@0:                     masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0:                     edgeweightlist = edgeweightlist(hitplanes);
tp@0:                     lastedgenumbers = lastedgenumbers(hitplanes);
tp@0:                     reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0:                     toedgecoords = toedgecoords(hitplanes,:);
tp@0: 	
tp@0:                     for ll = 1:jj
tp@0:                         eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0:                         if ll > kk
tp@0:                             eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0:                         end
tp@0:                     end
tp@0:                     ntocheckmore = length(masterlisttocheckmore);
tp@0:             
tp@0:                 end
tp@0:                 
tp@0:                 if SHOWTEXT >= 3
tp@0:                     disp(['         ',int2str(ntocheckmore),' of the special IS - N-edge - IR combinations survived the visibility test in refl plane ',int2str(kk)])
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:             
tp@0:             % Obstruction test of all the involved paths: R ->
tp@0:             % reflplane1 -> reflplane2 -> ... -> last edge
tp@0:             
tp@0:             if obstructtestneeded & ntocheckmore > 0
tp@0: 	
tp@0:                 for kk = 1:jj+1
tp@0:                     if ntocheckmore > 0
tp@0:                         if kk == 1
tp@0:                             fromcoords = R;
tp@0:                             startplanes = [];
tp@0:                         else
tp@0:                             eval(['fromcoords = reflpoints',JJ(kk-1,1:JJnumbofchars(kk-1)),';']) 
tp@0:                             startplanes = reflplanesexpand(:,kk-1);
tp@0:                         end
tp@0:                         if kk == jj+1,                    
tp@0:                             tocoords = toedgecoords;
tp@0:                             endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0:                         else
tp@0:                             eval(['tocoords = reflpoints',JJ(kk,1:JJnumbofchars(kk)),';'])    
tp@0:                             endplanes = reflplanesexpand(:,kk);    
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:                             masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0:                             edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0:                             lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0:                             reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0:                             toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0:                             for ll = 1:jj
tp@0:                                 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0:                                 eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0:                             end
tp@0:                             
tp@0:                         end
tp@0:                         ntocheckmore = length(masterlisttocheckmore);
tp@0:                         
tp@0:                     end
tp@0:                     
tp@0:                 end
tp@0:                 if SHOWTEXT >= 3
tp@0:                     disp(['         ',int2str(ntocheckmore),' of the special IS - N-edge - IR combinations survived the obstruction test'])
tp@0:                 end
tp@0:                 
tp@0:             end
tp@0:             
tp@0:             % Add the found special combinations to the outdata list
tp@0:         
tp@0:             edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,ii+1:ii+Ndifforder))-nplanes];
tp@0:             prespeclist = [prespeclist;[double(POTENTIALISES(masterlisttocheckmore,1:ii)) zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0:             postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0:             bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];    
tp@0:             eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(jj,1:JJnumbofchars(jj)),'];']);
tp@0:             % NB!! In the same way as earlier, we must a recursive reference
tp@0:             % method to find the image source of the last specular reflection.
tp@0:             ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterlisttocheckmore)));
tp@0:             for kk = 2:jj
tp@0:                 ivref = ORIGINSFROM(ivref);
tp@0:             end
tp@0:             validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0:             
tp@0:         end
tp@0:         
tp@0:         masterivlist = masterivlist(ivcompletelyOK);    
tp@0:         possibleedgecombs = possibleedgecombs(ivcompletelyOK,:);
tp@0:         possibleprespecs = possibleprespecs(ivcompletelyOK,:);
tp@0:         possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0:         possibleweights = possibleweights(ivcompletelyOK,:);
tp@0:             
tp@0:         nposs = length(ivcompletelyOK);
tp@0: 	
tp@0:         if SHOWTEXT >= 3
tp@0: 			disp(['         ',int2str(nposs),' IS + Edge (Nth-order) + IR segments (non-special combinations) survived the obstruction test'])
tp@0:         end
tp@0: 
tp@0:         % Add the found "standard" combinations to the outdata list
tp@0:     
tp@0:         edgedifflist = [edgedifflist;possibleedgecombs];
tp@0:         prespeclist = [prespeclist;  [possibleprespecs  zeros(nposs,specorder-Ndifforder-ii)]];
tp@0:         postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-Ndifforder-jj)]];
tp@0:         bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];    
tp@0:         validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0:         % NB!! In the same way as earlier, we must a recursive reference
tp@0:         % method to find the image source of the last specular reflection.
tp@0:         ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist)))));
tp@0:         for kk = 6:jj+Ndifforder
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: 
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: addvec = [];
tp@0: for ii = 1:Ndifforder
tp@0:     addvec = [addvec,zeros(length(edgedifflist),1),ones(length(edgedifflist),1)];
tp@0: end
tp@0: startandendpoints = [startandendpoints addvec];
tp@0: 
tp@0: %   #######################################################################
tp@0: %   #
tp@0: %   #   Construct a list guide, which will tell which rows have only
tp@0: %   #   dd, which rows have sdd etc
tp@0: %   #   Syntax:        dd,sdd,ssdd,sssdd,...,dds,ddss,ddsss,...
tp@0: %   #       (should also continue with sdds, sddss,...
tp@0: %   #   
tp@0: %   #######################################################################
tp@0: 
tp@0: [n1,n2] = size(prespeclist);
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 + Ndifforder - 1;
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,:);