tp@0: function edirfile = EDB1makeirs(edpathsfile,specorder,Rstart,...
tp@0:                         EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0:                         edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,...
tp@0:                         thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,...
tp@0:                         saveindividualdiffirs)
tp@0: % EDB1makeirs - Constructs impulse responses from a list of paths in a file.
tp@0: %
tp@0: % Input parameters:
tp@0: % 	edpathsfile         The name of the file containing the plane and edge data
tp@0: % 	specorder           The maximum order of reflections that is calculated.
tp@0: % 	Rstart              The reference distance that the time zero of the impulse
tp@0: % 	                    response refers to.
tp@0: % 	EDcalcmethod        'n' or 'v' or 'k', meaning:
tp@0: %                       'n': the new method by Svensson et al
tp@0: %                       'v': Vanderkooys method
tp@0: %                       'k': the Kirchhoff diffraction approximation
tp@0: % 	edgestartcoords     A matrix [nedges,3] of the startpoint coordinates
tp@0: % 	of
tp@0: % 	                    the edges.
tp@0: % 	edgeendcoords       A matrix [nedges,3] of the endpoint coordinates of
tp@0: % 	                    the edges.
tp@0: % 	edgenvecs           A matrix [nedges,3] of the normal vectors of
tp@0: % 	                    the reference plane of each edge.
tp@0: % 	edgelengthvec       A list [nedge,1] of the lenghts of each edge.
tp@0: %   planeeqs            A matrix, [nplanes,4], of all the plane equations.
tp@0: %   approxplanemidpoints   A matrix, [nplanes,3], of all plane midpoint
tp@0: %                       coordinates. It may be an approximate location of
tp@0: %                       the midpoint.
tp@0: % 	reflfactors         A list [nplanes,1] of the reflection factors of
tp@0: % 	                    each plane.
tp@0: % 	closwedangvec       A list [nedge,1] of the close-wedge angle of each
tp@0: % 	edge.
tp@0: % 	planesatedge        A matrix [nedge,2] of the two planes of each edge
tp@0: % 	                    with the reference plane first.
tp@0: % 	elemsize            A list [1,difforder] with the relative edge element
tp@0: % 	                    sizes that will be used for each order. The first
tp@0: % 	                    value, for diffraction order 1, is redundant and
tp@0: % 	                    not used. For higher orders, a value of, e.g., 0.5
tp@0: % 	                    means that edges will be subdivided into elements
tp@0: % 	                    such that the wavelength at half the sampling
tp@0: % 	                    frequency is 0.5 times the element size. A lower
tp@0: % 	                    value gives faster but less accurate results.
tp@0: % 	reftoshortlistE     A matrix, [nedges,nedges] with a pointer to the
tp@0: % 	                    short lists that contain edge-to-edge data.
tp@0: % 	re1sho              A short list, [nshort,1] of the cylindrical radii
tp@0: %                       of each edge startpoint relative to all other
tp@0: %                       edges.
tp@0: % 	re2sho              A short list, [nshort,1] of the cylindrical radii
tp@0: %                       of each edge endpoint relative to all other
tp@0: %                       edges.
tp@0: % 	thetae1sho          A short list, [nshort,1] of the theta angle
tp@0: %                       of each edge startpoint relative to all other
tp@0: %                       edges.
tp@0: % 	thetae2sho          A short list, [nshort,1] of the theta angle
tp@0: %                       of each edge endpoint relative to all other
tp@0: %                       edges.
tp@0: % 	ze1sho              A short list, [nshort,1] of the z-value
tp@0: %                       of each edge startpoint relative to all other
tp@0: %                       edges.
tp@0: % 	ze2sho              A short list, [nshort,1] of the z-value
tp@0: %                       of each edge endpoint relative to all other
tp@0: %                       edges.
tp@0: % 	edgeseespartialedge A matrix, [nedges,nedges], with edge-to-edge
tp@0: % 	                    visibility values. The values are 0 (invisible)
tp@0: %                       up to (2^(nedgesubs)-1)^2 (full visibility).
tp@0: %                       A pos. value indicates that the edge-to-edge path
tp@0: %                       runs along a plane; a neg. avlue indicates that it
tp@0: %                       does not run along a plane.
tp@0: %   edgeplaneperptoplane1   A matrix, [nplanes,nedges], with 0 or 1. A
tp@0: %                       value 1 indicates that one of the edge's two
tp@0: %                       defining planes is perpendicular to another plane.
tp@0: % 	desiredirfile       The file name that will be given to the
tp@0: % 	                    output file. 
tp@0: %   guiderowstouse      (optional) A list of values, indicating which rows in the
tp@0: %                       mainlistguide and mainlistguidepattern that should
tp@0: %                       be used. This way it is possible to select only
tp@0: %                       diffraction, for instance. If this list is not
tp@0: %                       specified, all components will be used.
tp@0: %   includedirectsound  (optional) 0 or 1, indicating whether the direct
tp@0: %                       sound should be included or not. Default: 1 
tp@0: %   saveindividualdiffirs (optional) Vector of length 2 with values 0 or 1
tp@0: %						Pos 1 indicating whether
tp@0: %                       	0 - all orders of diffraction IR:s will be summed
tp@0: %                       	in a single vector 'irdiff', or
tp@0: %                       	1 - each order of diffraction IR:s will be placed
tp@0: %                       	in individual columns in a matrix 'irdiff'.
tp@0: %						Pos 2 indicating whether
tp@0: %							0 - only the total diffraction ir will be saved in the file.
tp@0: %							1 - all individual diffraction irs will be saved in a large 
tp@0: %							matrix alldiffirs.
tp@0: %                       Default: [0 0].
tp@0: %   FSAMP, CAIR, SHOWTEXT   Global parameters.
tp@0: %
tp@0: % Output parameters:
tp@0: %   edirfile            The filename used for the output file that contains
tp@0: %                       the impulse responses.
tp@0: %
tp@0: % Data in the output file:
tp@0: %   irdirect            The IR containing the direct sound, size [nirlength,1].
tp@0: %   irgeom              The IR containing the specular reflections, size [nirlength,1].
tp@0: %   irdiff              The IR containing the diffracted components, size [nirlength,1].
tp@0: %   irtot               The IR containing the sum of all components, size [nirlength,1].
tp@0: %	alldiffirs			(optional) Matrix with all individual first-order diffraction IRs
tp@0: %						on individual lines.
tp@0: %	alldiffirsextradata	(optional) Vector with combination number (form reflpaths) that
tp@0: %						matches alldiffirs.
tp@0: %   FSAMP Rstart CAIR   Same as the input parameters.
tp@0: %
tp@0: % Uses functions  EDB1irfromslotlist EDB1calcdist EDB1coordtrans1 EDB1coordtrans2
tp@0: % 				  EDB1wedge1st_int EDB1wedge2nd
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) 20061208
tp@0: %
tp@0: % edirfile = EDB1makeirs(edpathsfile,specorder,...
tp@0: %            Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0: %            edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,...
tp@0: %            thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,...
tp@0: %            saveindividualdiffirs);
tp@0: 
tp@0: global SHOWTEXT FSAMP CAIR BIGEDGESTEPMATRIX
tp@0: 
tp@0: global IRDIFFVEC
tp@0: 
tp@0: eval(['load ',edpathsfile])
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: % We look for the optional vector multfactors in the edpathsfile.
tp@0: % If it was there, then those values will be used to switch off, or boost
tp@0: % any reflection path.
tp@0: % Those values will only be used for diffraction paths.
tp@0: 
tp@0: if exist('multfactors') == 0
tp@0:     multfactors = ones(size(reflpaths,1),1);    
tp@0: else
tp@0:     if size(multfactors,1) ~= size(reflpaths,1)
tp@0:         error(['The edpaths file contained a vector multfactors which does not have the same length as reflpaths.'])    
tp@0:     end
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------------------------
tp@0: 
tp@0: edirfile = desiredirfile;
tp@0: 
tp@0: Sdirection = [1 0 0];
tp@0: maxnedgesorplanes = max(max(reflpaths));
tp@0: if ~isempty(maxnedgesorplanes)
tp@0:     multfac = 10^ceil(log10(double(maxnedgesorplanes)+1));
tp@0: else
tp@0:     multfac = 0;
tp@0: end
tp@0: 
tp@0: if nargin < 27
tp@0:     saveindividualdiffirs = [0 0];
tp@0: end
tp@0: if nargin < 26
tp@0:     includedirectsound = 1;
tp@0: end
tp@0: if nargin < 25
tp@0:     guiderowstouse = [];
tp@0: end
tp@0: if isempty(guiderowstouse)
tp@0:     usesubset = 0;    
tp@0:     if SHOWTEXT > 2
tp@0:         disp('Using all components')
tp@0:     end
tp@0: else
tp@0:     usesubset = 1;        
tp@0:     if SHOWTEXT > 2
tp@0:         disp('Using only some components')
tp@0:     end
tp@0: end
tp@0: 
tp@0: nyvec = pi./(2*pi - closwedangvec);    
tp@0: onesvec1 = ones(1,3);
tp@0: 
tp@0: souspecboost = 1;
tp@0: if ~isempty(Sinsideplanenumber)
tp@0:     if reflfactors(Sinsideplanenumber(1)) ~= 0
tp@0:         souspecboost = 2;
tp@0:     end
tp@0: end
tp@0: 
tp@0: recspecboost = 1;
tp@0: if ~isempty(Rinsideplanenumber)
tp@0:     if reflfactors(Rinsideplanenumber(1)) ~= 0
tp@0:         recspecboost = 2;
tp@0:     end
tp@0: end
tp@0: 
tp@0: if isempty(re1sho)
tp@0:     userwantsdiff2 = 0;
tp@0: else
tp@0:     userwantsdiff2 = 1;    
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: 
tp@0: if exist('mainlistguide') ~= 1
tp@0:     mainlistguide = [];    
tp@0: else
tp@0:     mainlistguide = double(mainlistguide);    
tp@0: end
tp@0: 
tp@0: [ncomponents,ncols] = size(reflpaths);
tp@0: [nrowsinlist,slask] = size(mainlistguide);
tp@0: if ncols > 1
tp@0:     difforderinlist = sum(mainlistguidepattern.'=='d').';
tp@0: else
tp@0:     difforderinlist = (mainlistguidepattern=='d');    
tp@0: end
tp@0: 
tp@0: lastNdiffrow = zeros(1,specorder);
tp@0: for ii = 1:specorder
tp@0:     iv = find(difforderinlist <= ii );
tp@0:     if ~isempty(iv)
tp@0:         lastNdiffrow(ii) = iv(end);
tp@0:     end
tp@0: end
tp@0: 
tp@0: nrefltogetherwdiff = sum(mainlistguidepattern.'=='d').'.*(  sum(mainlistguidepattern.'=='s').');
tp@0: 
tp@0: if SHOWTEXT >= 3
tp@0: 	disp(['   Constructing IR. ',int2str(ncomponents),' components.'])
tp@0: end
tp@0: 
tp@0: irdirect = [0 0].';
tp@0: irgeom = [0 0].';
tp@0: irdiff = [0 0].';
tp@0: irtot = [0 0].';
tp@0: 
tp@0: Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod'];
tp@0: 
tp@0: if ncomponents == 0
tp@0:     eval(['save ',edirfile,Varlist])
tp@0:     return
tp@0: end
tp@0: 
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %
tp@0: %  Diffraction once, possibly with pre- and post-specular reflections
tp@0: %
tp@0: %##############################################################################
tp@0: 
tp@0: directsoundonboundary = 0;
tp@0: specreflonboundary = zeros(size(planeeqs,1),1);
tp@0: 
tp@0: if firstdiffrow ~= 0
tp@0:     
tp@0:     % First we remove the 'd' for all the rows in the mainlistguidepattern that
tp@0:     % we do not want to use.
tp@0:     
tp@0:     if usesubset == 1
tp@0:         singdiffrows = find(difforderinlist==1);
tp@0:         rowstoclear = setdiff(singdiffrows,guiderowstouse);
tp@0:         if ~isempty(rowstoclear)
tp@0:             mainlistguidepattern(rowstoclear,:) = mainlistguidepattern(rowstoclear,:)*0;    
tp@0:         end
tp@0:     end
tp@0: 
tp@0:     % Here we can use the single-diffraction calculation for all
tp@0:     % combinations. The source should either be the original source or
tp@0:     % an image source. The receiver should either be the original
tp@0:     % receiver or an image receiver.
tp@0:     
tp@0:     % First we create lists that specify whether there are specular
tp@0:     % reflections before the edge diffraction, and after the edge
tp@0:     % diffraction. Also, a list which gives the column where the edge
tp@0:     % diffraction is present can be used to efficiently extract the
tp@0:     % edge number. 
tp@0:     % NB! The lists prespeclist, postspeclist and singdiffcol all have
tp@0:     % the (short) length of the mainlistguide. 
tp@0:     % The data in these lists will then be used to create the long
tp@0:     % lists prespecpresent and postspecpresent which simply are 0 or 1.
tp@0:     % A matrix called edgemask will have the same number of columns as
tp@0:     % the reflpaths matrix. For each row, there will be a 1 in the
tp@0:     % column where the edge diffraction is so we can efficiently find
tp@0:     % the edge number.
tp@0: 
tp@0:     multmatrix = [1:specorder];
tp@0:     if size(mainlistguidepattern,2) > specorder
tp@0:         multmatrix = [1:size(mainlistguidepattern,2)];
tp@0:     end
tp@0:     multmatrix = multmatrix( ones(nrowsinlist,1),:  );
tp@0:     
tp@0:     if ncols > 1
tp@0:         singdiffcol = sum( (multmatrix.*(mainlistguidepattern=='d')).' ).';
tp@0:         singdiffcol = singdiffcol.*(sum(mainlistguidepattern.'=='d').'<=1);
tp@0:         nrefl = sum( (mainlistguidepattern=='s' | mainlistguidepattern=='d').' ).';
tp@0:     else
tp@0:         singdiffcol = mainlistguidepattern=='d';
tp@0:         nrefl = singdiffcol | mainlistguidepattern=='s';
tp@0:     end
tp@0:     
tp@0:     prespeclist = singdiffcol-1;
tp@0:     prespeclist = prespeclist.*(prespeclist>0);
tp@0:     postspeclist = (nrefl-singdiffcol).*(singdiffcol~=0);
tp@0:     
tp@0:     prespecpresent = zeros(ncomponents,1);
tp@0:     postspecpresent = zeros(ncomponents,1);
tp@0:     edgemask = zeros(ncomponents,specorder);
tp@0: 
tp@0:     diffrowsthatareOK = [firstdiffrow:lastNdiffrow(1)];
tp@0:     if usesubset == 1
tp@0:         diffrowsthatareOK = intersect(diffrowsthatareOK,guiderowstouse);
tp@0:     end
tp@0:     for ii = diffrowsthatareOK,    %firstdiffrow:lastNdiffrow(1)
tp@0:         iv = [(mainlistguide(ii,2)):(mainlistguide(ii,3))];
tp@0:         onesvec2 = ones(length(iv),1);
tp@0:         edgemask(iv,singdiffcol(ii)) = onesvec2;            
tp@0:     	prespecpresent(iv) = (prespeclist(ii)>0)*onesvec2;
tp@0:         postspecpresent(iv) = (postspeclist(ii)>0)*onesvec2;
tp@0:     end
tp@0:     prespecpresent = prespecpresent(:,onesvec1);
tp@0:     postspecpresent = postspecpresent(:,onesvec1);
tp@0: 
tp@0:     if ncols > 1
tp@0:         edgenumb = sum( (edgemask.*double(reflpaths(:,1:specorder))).').';
tp@0:     else
tp@0:         edgenumb = edgemask.*double(reflpaths);    
tp@0:     end
tp@0: 
tp@0:     nedgeirs = size(edgenumb,1);
tp@0:     nnonzeroedgeirs = sum(edgenumb(:,1)>0);
tp@0: 
tp@0:     ircounter = 0;
tp@0:     edgelist = unique(edgenumb);
tp@0:     for ii = 1: length(edgelist)
tp@0:         edge = edgelist(ii);
tp@0:         if edge~= 0
tp@0:            iv = find(edgenumb==edge);
tp@0:            ncombs = length(iv);
tp@0:            
tp@0:            if ncombs > 0 & any(multfactors(iv))
tp@0:                if SHOWTEXT >= 4
tp@0:                     disp(['   Edge ',int2str(edge),': ',int2str(ncombs),' combinations'])    
tp@0:                end               
tp@0:                  
tp@0:                 onesvec3 = ones(ncombs,1);
tp@0: 	
tp@0:                 IS = full(specextradata(iv,1:3));
tp@0:                 IR = full(specextradata(iv,4:6));                
tp@0:             
tp@0:                 edgestart = full(edgeextradata(iv,1));
tp@0:                 edgeend   = full(edgeextradata(iv,2));
tp@0:                 
tp@0:                 % Calculate the cyl coordinates of all IS
tp@0:                 % 
tp@0:                 % Calculate the cyl coord of all IR
tp@0: 	
tp@0:                 edgecoords = [edgestartcoords(edge,:);edgeendcoords(edge,:)];
tp@0:             	[rs,thetas,zs,rr,thetar,zr] = EDB1coordtrans2(IS,IR,edgecoords,edgenvecs(edge,:));
tp@0: 	
tp@0: 				bc = reflfactors(planesatedge(edge,:)).';
tp@0: 				if prod(double(bc==[1 1])) ~= 1
tp@0: 					disp('   Non-rigid wedge')
tp@0: 				end
tp@0: 	
tp@0:                 for jj = 1:ncombs
tp@0:                     if multfactors(iv(jj)) > 0
tp@0:                         if EDcalcmethod(1) == 'n'
tp@0: 
tp@0:                            [irnew,slask,singularterm] = EDB1wedge1st_int(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),...
tp@0:                              edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc);                  
tp@0:                          
tp@0:                               if any(singularterm)
tp@0:                                   if singularterm(2) | singularterm(3)
tp@0:                                       directsoundonboundary = 1;
tp@0:                                 elseif singularterm(1)
tp@0:                                     if specorder == 1
tp@0:                                        specreflonboundary(planesatedge(edge,2)) = 1;
tp@0:                                    else
tp@0:                                        disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!');    
tp@0:                                    end
tp@0:                                 elseif singularterm(4)
tp@0:                                       if specorder == 1
tp@0:                                            specreflonboundary(planesatedge(edge,1)) = 1;
tp@0:                                        else
tp@0:                                            disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!');    
tp@0:                                        end
tp@0:                                 end
tp@0:                               end
tp@0: 		
tp@0:                         else    % ...   if EDcalcmethod(1) == 'n'
tp@0:                            
tp@0:                               [irnew,slask,singularterm] = EDB1wedge1stcombo(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),...
tp@0:                                   edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc);                      
tp@0:                               
tp@0:                         end      % ...   if EDcalcmethod(1) == 'n'
tp@0:                     
tp@0:                         % Decide whether the IR should be boosted or not
tp@0:                         %
tp@0:                         % The factors souspecboost and recspecboost have
tp@0:                         % values other than 1 if the source or receiver is directly
tp@0:                         % at a plane, for the boosting of the direct sound, and
tp@0:                         % specular reflections.
tp@0:                         %
tp@0:                         % This boost factor should be used also for edge
tp@0:                         % diffraction if:
tp@0:                         %   1.  There are some specular reflections between the 
tp@0:                         %       source and the edge
tp@0:                         %   or
tp@0:                         %   2.  There are no specular reflections between the
tp@0:                         %       source and the edge, but the source/receiver is at
tp@0:                         %       a plane which does not belong to the edge.
tp@0:                         
tp@0:                         if souspecboost ~= 1 | recspecboost ~= 1
tp@0:                         
tp@0:                             boostfactor = 1;
tp@0:                             if prespecpresent(iv(jj)) == 1
tp@0:                                 boostfactor = souspecboost;    
tp@0:                             else
tp@0:                                 if ~isempty(Sinsideplanenumber)
tp@0:                                     if Sinsideplanenumber(1)~=planesatedge(edge,1) & Sinsideplanenumber(1)~=planesatedge(edge,2), 
tp@0:                                         boostfactor = souspecboost;
tp@0:                                     end
tp@0:                                 end
tp@0:                             end
tp@0:                             if postspecpresent(iv(jj)) == 1
tp@0:                                 boostfactor = boostfactor*recspecboost;    
tp@0:                             else
tp@0:                                 if ~isempty(Rinsideplanenumber)
tp@0:                                     if Rinsideplanenumber(1)~=planesatedge(edge,1) & Rinsideplanenumber(1)~=planesatedge(edge,2), 
tp@0:                                         boostfactor = boostfactor*recspecboost;
tp@0:                                     end
tp@0:                                 end
tp@0:                             end
tp@0:                             if boostfactor ~= 1
tp@0:                                 irnew = irnew*boostfactor;    
tp@0:                             end
tp@0:                         
tp@0:                         end    % ....    if souspecboost ~= 1 | recspecboost ~= 1
tp@0: 
tp@0:                         ndiff = length(irdiff);
tp@0:                         nnew = length(irnew);
tp@0:                         
tp@0:                         ircounter = ircounter + 1;                        
tp@0:                         
tp@0:                         if nnew > ndiff
tp@0:                            irdiff = [irdiff;zeros(nnew-ndiff,1)];
tp@0:                         end
tp@0:                         irdiff(1:nnew) = irdiff(1:nnew) + irnew*double(multfactors(iv(jj)));
tp@0: 
tp@0: 						if saveindividualdiffirs(2) == 1
tp@0: 							if exist('alldiffirs','var') == 0
tp@0: 								alldiffirs = sparse(zeros(nnonzeroedgeirs,nnew));
tp@0: 								alldiffirs(1,:) = irnew*double(multfactors(iv(jj))).';
tp@0: 								alldiffirsextradata = zeros(nnonzeroedgeirs,1);
tp@0: 								alldiffirsextradata(1) = iv(jj);
tp@0: 							else
tp@0: 								if nnew > ndiff
tp@0: 									alldiffirs = [alldiffirs zeros(nnonzeroedgeirs,nnew-ndiff)];
tp@0: 								end
tp@0: 								alldiffirs(ircounter,1:nnew) = irnew*double(multfactors(iv(jj))).';							
tp@0: 								alldiffirsextradata(ircounter) = iv(jj); 
tp@0: 							end						
tp@0: 						
tp@0: 						end
tp@0: 
tp@0: 
tp@0:                         if SHOWTEXT >= 7
tp@0:                             figure(1)
tp@0:                             plot(irnew)
tp@0:                             figure(2)
tp@0:                             plot(irdiff)
tp@0:                             save ~/Documents/Temp/irdiff2.mat
tp@0:                             pause
tp@0:                         end  
tp@0: 
tp@0:                     end   %  ...  if multfactors(iv) > 0
tp@0:                     
tp@0:                 end  % .....                for jj = 1:ncombs
tp@0: 
tp@0:             end   %....    if ncombs > 0 & any(multfactors(iv))
tp@0:             
tp@0:         end    % ....         if edge~= 0
tp@0:         
tp@0:     end  % ....     for ii = 1: length(edgelist)
tp@0:     
tp@0: end   % ...if firstdiffrow ~= 0
tp@0: 
tp@0: nspecreflonboundary = sum(specreflonboundary);
tp@0: 
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %
tp@0: %      The direct sound
tp@0: %
tp@0: %##############################################################################
tp@0: 
tp@0: if usesubset == 1 & directsoundrow == 1
tp@0:     if any(guiderowstouse == 1) == 0
tp@0:         directsoundrow = 0;    
tp@0:     end
tp@0: end
tp@0: 
tp@0: if directsoundrow == 1 & includedirectsound == 1
tp@0:     dist = norm(R-S);
tp@0:     slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1;
tp@0:     amp = souspecboost*recspecboost./dist;
tp@0:     if directsoundonboundary
tp@0:         amp = amp/2;    
tp@0:         if SHOWTEXT >= 4
tp@0:             disp('HALVING DIRECT SOUND')
tp@0:         end
tp@0:     end
tp@0:     slotnumber = floor(slotnumberfrac);
tp@0:     if any(slotnumber<1)
tp@0:         error('ERROR: Rstart is set to too low a value!')
tp@0:     end
tp@0:     slotnumberfrac = slotnumberfrac - slotnumber;
tp@0:     irdirect = zeros( slotnumber+2,1 );
tp@0: 
tp@0: 	irdirect(slotnumber)   = amp.*(1-slotnumberfrac) ;
tp@0: 	irdirect(slotnumber+1) = amp.*slotnumberfrac;	
tp@0: 
tp@0: 	if ncomponents == 1
tp@0:         nnew = length(irdirect);
tp@0:         if nnew > 2
tp@0:           irdiff = [irdiff;zeros(nnew-2,1)];
tp@0:           irtot =  [irtot;zeros(nnew-2,1)];
tp@0:           irgeom = [irgeom;zeros(nnew-2,1)];
tp@0:         end
tp@0:         irtot(1:nnew) = irtot(1:nnew) + irdirect;
tp@0: 
tp@0:         irtot = sparse(irtot);
tp@0:         irgeom = sparse(irgeom);
tp@0:         irdirect = sparse(irdirect);
tp@0:         irdiff = sparse(irdiff);
tp@0:         
tp@0:         eval(['save ',edirfile,Varlist])
tp@0:         return
tp@0: 	end    
tp@0: else
tp@0:     if directsoundonboundary == 1
tp@0:         
tp@0:         dist = norm(R-S);
tp@0:         slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1;
tp@0:         amp = souspecboost*recspecboost./dist;
tp@0:         if directsoundonboundary
tp@0:             amp = amp/2;    
tp@0:             if SHOWTEXT >= 4
tp@0:                 disp('HALVING DIRECT SOUND')
tp@0:             end
tp@0:         end
tp@0:         slotnumber = floor(slotnumberfrac);
tp@0:         if any(slotnumber<1)
tp@0:             error('ERROR: Rstart is set to too low a value!')
tp@0:         end
tp@0:         slotnumberfrac = slotnumberfrac - slotnumber;
tp@0:         irdirect = zeros( slotnumber+2,1 );
tp@0: 	
tp@0: 		irdirect(slotnumber)   = amp.*(1-slotnumberfrac) ;
tp@0: 		irdirect(slotnumber+1) = amp.*slotnumberfrac;	
tp@0:         
tp@0:     end
tp@0:     if SHOWTEXT >= 4
tp@0:         disp(['      No direct sound'])    
tp@0:     end
tp@0: end
tp@0: 
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %
tp@0: %  All-specular			s, ss, sss, etc
tp@0: %
tp@0: %##############################################################################
tp@0: 
tp@0: if allspecrows(1) ~= 0
tp@0:     if usesubset == 1
tp@0:         specrowsthatareOK = intersect(allspecrows,guiderowstouse);
tp@0:     else
tp@0:         specrowsthatareOK = [allspecrows(1):allspecrows(2)];
tp@0:     end
tp@0:     if ~isempty(specrowsthatareOK)
tp@0:         if usesubset == 1
tp@0:             temp = mainlistguide(specrowsthatareOK,2:3);
tp@0:             ivspec = [double(mainlistguide(specrowsthatareOK(1),2)):double(mainlistguide(specrowsthatareOK(1),3))];
tp@0:             for ii = 2:size(temp,1);
tp@0:                 ivspec = [ivspec  [double(mainlistguide(specrowsthatareOK(ii),2)):double(mainlistguide(specrowsthatareOK(ii),3))]];
tp@0:             end
tp@0:         else
tp@0:             ivspec = [mainlistguide(allspecrows(1),2):mainlistguide(allspecrows(2),3)];        
tp@0:         end
tp@0:         
tp@0:         if nspecreflonboundary > 0
tp@0:             listofspecreflonboundary = find(specreflonboundary);
tp@0:             specrefltoadd = setdiff(listofspecreflonboundary,ivspec)
tp@0:             if ~isempty(specrefltoadd)
tp@0:                 error(['ERROR: We need to add some code here, to reinsert pruned spec. refl.'])    
tp@0:             end
tp@0:         end
tp@0:         
tp@0:         if SHOWTEXT >= 4
tp@0:             disp(['      ',int2str(length(ivspec)),' specular reflections'])    
tp@0:         end
tp@0: 	   
tp@0:         dist = EDB1calcdist(full(specextradata(ivspec,1:3)),R);
tp@0:     
tp@0: 		if nspecreflonboundary > 0
tp@0: 			specamp = ones(size(souspecboost));
tp@0: 			amplitudeshouldbehalf = ismember(reflpaths(ivspec,1),listofspecreflonboundary);
tp@0: 			specamp = specamp - amplitudeshouldbehalf*0.5;
tp@0: 			irnew = specamp*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);		                
tp@0: 		else
tp@0: 			irnew = souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);		
tp@0: 		end
tp@0: 		ngeom = length(irgeom);
tp@0: 		nnew = length(irnew);
tp@0: 		if nnew > ngeom
tp@0: 		   irgeom = [irgeom;zeros(nnew-ngeom,1)];
tp@0: 		end
tp@0: 		irgeom(1:nnew) = irgeom(1:nnew) + irnew;
tp@0: 		
tp@0:     else
tp@0:         if nspecreflonboundary > 0
tp@0:             listofspecreflonboundary = find(specreflonboundary);
tp@0:             error(['ERROR: We need to add some code here, to reinsert pruned spec. refl., pos. 2'])    
tp@0:         end        
tp@0:     end    
tp@0: else
tp@0:     if nspecreflonboundary > 0
tp@0:         listofspecreflonboundary = find(specreflonboundary);
tp@0:         
tp@0:         [xis] = EDB1findis(S,listofspecreflonboundary,planeeqs,1,[1 1 1]);
tp@0:         
tp@0:         disp(['WARNING: We have some new code here, to reinsert pruned spec. refl., pos. 3'])
tp@0:         dist = EDB1calcdist(xis,R);
tp@0:         irnew = 0.5*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);		                
tp@0:         ngeom = length(irgeom);
tp@0:         nnew = length(irnew);
tp@0:         if nnew > ngeom
tp@0:            irgeom = [irgeom;zeros(nnew-ngeom,1)];
tp@0:         end
tp@0:         irgeom(1:nnew) = irgeom(1:nnew) + irnew;
tp@0:         
tp@0:     end        
tp@0: end
tp@0: 
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %##############################################################################
tp@0: %
tp@0: %  Multiple diffraction, possibly with pre- and post-specular reflections
tp@0: %
tp@0: %##############################################################################
tp@0: 
tp@0: if userwantsdiff2 == 1
tp@0:     
tp@0:     JJ = setstr(32*ones(specorder,1));
tp@0:     for jj=1:specorder
tp@0:         jjstr = int2str(jj);
tp@0:         JJ(jj,1:length(jjstr)) = jjstr;
tp@0:     end
tp@0:     [n1,n2] = size(JJ);
tp@0: 
tp@0:     
tp@0: 	for Ndifforder = 2:specorder
tp@0:         if SHOWTEXT >= 3
tp@0:             disp(['   Diffraction order ',int2str(Ndifforder)])    
tp@0:         end
tp@0: 	
tp@0:         if any(difforderinlist==Ndifforder) & elemsize(Ndifforder) > 0
tp@0: 	
tp@0:             % Calculate some general parameters that are shared for all
tp@0:             % N-diffraction calculations
tp@0:             
tp@0:             divmin = CAIR/(FSAMP*elemsize(Ndifforder));
tp@0:             ndivvec = ceil(abs( edgelengthvec.' )/divmin);
tp@0:             dzvec = (edgelengthvec.')./ndivvec;
tp@0: 	
tp@0:             ncylrows = 4*(Ndifforder-1);
tp@0:             
tp@0:             % Here we can use the double-diffraction calculation for all
tp@0:             % combinations. The source should either be the original source or
tp@0:             % an image source. The receiver should either be the original
tp@0:             % receiver or an image receiver.
tp@0:             
tp@0:             noffset = lastNdiffrow(Ndifforder-1);
tp@0:             ivNdiff = [noffset+1:lastNdiffrow(Ndifforder)];
tp@0:             ndiffcombs = length(ivNdiff);
tp@0:             zerosvec1 = zeros(ndiffcombs,1);
tp@0:             
tp@0:             ndoubrows = length(ivNdiff);        
tp@0:             previousrow = lastNdiffrow(Ndifforder-1);
tp@0:             if previousrow > 0
tp@0:                 noffsetlonglist = mainlistguide(previousrow,3);
tp@0:             else
tp@0:                 noffsetlonglist = 0;    
tp@0:             end
tp@0:             nremaining = ncomponents - (mainlistguide(lastNdiffrow(Ndifforder),3));
tp@0:             nlonglist = ncomponents-noffsetlonglist-nremaining;
tp@0:             ivlonglist = [mainlistguide(ivNdiff(1),2):mainlistguide(ivNdiff(end),3)].';
tp@0: 	
tp@0:             diffcols = zerosvec1(:,ones(1,Ndifforder));
tp@0:             for ii = ivNdiff(1):ivNdiff(end)
tp@0:                 diffcols(ii-ivNdiff(1)+1,:) = find(mainlistguidepattern(ii,:)=='d');
tp@0:             end
tp@0:             
tp@0:             nreflorder = sum( (mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='d'|mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='s').').';
tp@0:             nprespecs = diffcols(:,1)-1;
tp@0:             nmidspecs = diffcols(:,2)-diffcols(:,1)-1;
tp@0:             npostspecs = nreflorder-diffcols(:,Ndifforder);
tp@0: 	
tp@0:             if ndiffcombs > 1
tp@0:                 ivreftolonglist = ivlonglist(:,ones(1,ndiffcombs));
tp@0:                 comppattern = mainlistguide(ivNdiff,2).';        
tp@0:                 comppattern = comppattern(ones(nlonglist,1),:);
tp@0:             
tp@0:                 rowinpatternlist = sum( (ivreftolonglist>=comppattern).' ).';
tp@0:             else
tp@0:                 rowinpatternlist = ones(size(ivlonglist));    
tp@0:             end
tp@0:             
tp@0:             % Construct a long matrix with the edge numbers extracted from the
tp@0:             % reflpaths matrix.
tp@0: 	
tp@0:             longnmidspec = zeros(nlonglist,1);
tp@0:             iv = find(nmidspecs>0);
tp@0:             for ii = 1:length(iv)
tp@0:                 ivreftolonglist = [mainlistguide(ivNdiff(iv(ii)),2):mainlistguide(ivNdiff(iv(ii)),3)] - noffsetlonglist;
tp@0:                 longnmidspec(ivreftolonglist) = nmidspecs(iv(ii));
tp@0:             end
tp@0:             
tp@0:             %------------------------------------------------------------------
tp@0:             % edgepattern will be a matrix which, for each reflection combination, contains
tp@0:             % the 2,3,4,... edge numbers that are involved in each path
tp@0:             % regardless of there are specular reflections before, in the
tp@0:             % middle, or after.
tp@0:             
tp@0:             edgepattern = zeros(nlonglist,Ndifforder);
tp@0:             countvec = [1:nlonglist].';
tp@0:             reflpathscut = reflpaths(ivlonglist,:);            
tp@0:             for ii = 1:Ndifforder
tp@0:                 ivrefcol = countvec + (diffcols(rowinpatternlist,ii)-1)*nlonglist;
tp@0:                 edgepattern(:,ii) = reflpathscut(ivrefcol);
tp@0:             end
tp@0: 	
tp@0:             %------------------------------------------------------------------
tp@0:             % midspecpattern will be a matrix which, for each reflection combination, contains
tp@0:             % the 1,2,3,4,... specular reflection numbers (i.e., plane numbers) that are involved
tp@0:             % in each path between diffraction 1 and diffraction 2.
tp@0:             % First, the reflection component immediately after diffraction 1
tp@0:             % is selected for every reflection. Then there is a masking which
tp@0:             % zeroes the midspecpattern value for all the combinations that
tp@0:             % don't have any midspec.
tp@0:             % NB! For combinations like ssssdd, the first step will point to a
tp@0:             % column which is outside the matrix. This is fixed by maximizing
tp@0:             % the column number.
tp@0:             
tp@0:             midspecpattern = zeros(nlonglist,max([specorder-Ndifforder 1]));
tp@0:             maxpointvalue = nlonglist*max([specorder-Ndifforder 1]);
tp@0:             for ii = 1:specorder-Ndifforder
tp@0:                 ivrefcol = countvec + (diffcols(rowinpatternlist,1)+(ii-1))*nlonglist;
tp@0:                 ivrefcol = mod(ivrefcol-1,maxpointvalue)+1;
tp@0:                 midspecpattern(:,ii) = double(reflpathscut(ivrefcol)).*(longnmidspec>=ii);
tp@0:             end
tp@0:             diff1col = diffcols(rowinpatternlist,1);
tp@0: 	
tp@0:             % Many combinations include the same edge-pair/N-let, so we extract the
tp@0:             % unique edge pairs/N-lets, and go through them in the ii-loop
tp@0:             
tp@0:             [edgeshortlist,ivec,jvec] = unique(edgepattern,'rows');
tp@0:             lastndivcomb = zeros(1,Ndifforder);
tp@0:             
tp@0:             for ii = 1: length(ivec)
tp@0:                 
tp@0:                 if SHOWTEXT >= 3
tp@0:                     if round(ii/ceil(length(ivec)/10))*ceil(length(ivec)/10) == ii
tp@0:                         disp(['      Combination no. ',int2str(ii),' of ',int2str(length(ivec))]) 
tp@0:                     end
tp@0:                 end
tp@0:                 
tp@0:                 iv = find(jvec==ii);
tp@0:                 ncombs = length(iv);
tp@0:                 onesvec3 = ones(ncombs,1);
tp@0: 	
tp@0:                 if SHOWTEXT >= 4
tp@0:                     numvec = int2str(edgeshortlist(ii,1));
tp@0:                     for ll = 2:Ndifforder
tp@0:                         numvec = [numvec,' ',  int2str(edgeshortlist(ii,ll))];
tp@0:                     end
tp@0:                     disp(['   Edges ',numvec])
tp@0:                         
tp@0:                 end
tp@0:                 
tp@0:                 if Ndifforder >= 2 & any(multfactors(ivlonglist(iv)))
tp@0: 	
tp@0:                     newndivvec = ndivvec(edgeshortlist(ii,:));
tp@0:                     
tp@0:                     if any(lastndivcomb~=newndivvec)
tp@0:                         if Ndifforder == 2
tp@0:                             ivmatrix = EDB1creindexmatrix(newndivvec);
tp@0:                         else
tp@0:                             ivmatrix = EDB1creindexmatrix(newndivvec(2:end));                        
tp@0:                         end
tp@0:                         [nedgeelcombs,slask] = size(ivmatrix);
tp@0:                         if Ndifforder == 2
tp@0:                             BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),:);
tp@0:                         else
tp@0:                             BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),2:end);
tp@0:                         end
tp@0:                         clear ivmatrix
tp@0:                         lastndivcomb = newndivvec;
tp@0: 	
tp@0:                     end
tp@0:                 end            
tp@0: 
tp@0:                 pathalongplane = (edgeseespartialedge(edgeshortlist(ii,2),edgeshortlist(ii,1))>0);
tp@0:                 for jj = 3:Ndifforder
tp@0:                     pathalongplane = [pathalongplane,(edgeseespartialedge(edgeshortlist(ii,jj),edgeshortlist(ii,jj-1))>0)];                    
tp@0:                 end
tp@0:                 
tp@0:                 if any(multfactors(ivlonglist(iv)))
tp@0: 	
tp@0:                     IS = full(specextradata(ivlonglist(iv),1:3));
tp@0:                     IR = full(specextradata(ivlonglist(iv),4:6));
tp@0:                     
tp@0:                     firstedgestart = full(edgeextradata(ivlonglist(iv),1));
tp@0:                     firstedgeend   = full(edgeextradata(ivlonglist(iv),2));
tp@0: 		
tp@0:                     lastedgestart = full(edgeextradata(ivlonglist(iv),3));
tp@0:                     lastedgeend   = full(edgeextradata(ivlonglist(iv),4));
tp@0: 		
tp@0:                     % Calculate the cyl coordinates of all IS/S and IR/R
tp@0: 		
tp@0:                     cylS = zeros(ncombs,3);
tp@0:                     edgecoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)];
tp@0:                 	[cylS(:,1),cylS(:,2),cylS(:,3)] = EDB1coordtrans1(IS,edgecoords,edgenvecs(edgeshortlist(ii,1),:));               
tp@0: 		
tp@0:                     cylR = zeros(ncombs,3);
tp@0:                     edgecoords = [edgestartcoords(edgeshortlist(ii,Ndifforder),:);edgeendcoords(edgeshortlist(ii,Ndifforder),:)];
tp@0:                 	[cylR(:,1),cylR(:,2),cylR(:,3)] = EDB1coordtrans1(IR,edgecoords,edgenvecs(edgeshortlist(ii,Ndifforder),:));     
tp@0:                         
tp@0:                     bc = ones(Ndifforder,2);    % Check real reflfactors!!!    
tp@0:                     bc = reshape(bc.',2*Ndifforder,1);
tp@0: 		
tp@0:                     % Pick out the edge-to-edge coordinates for this specific
tp@0:                     % edge pair/N-let.                        
tp@0:                                                 
tp@0:                       if ~isempty(reftoshortlistE),                
tp@0:                             if Ndifforder >= 2
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,1));
tp@0:                                 cylE2_r1 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,1),edgeshortlist(ii,2));
tp@0:                                 cylE1_r2 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 3
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,2));
tp@0:                                 cylE3_r2 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,3));
tp@0:                                 cylE2_r3 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 4
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,3));
tp@0:                                 cylE4_r3 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,4));
tp@0:                                 cylE3_r4 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 5
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,4));
tp@0:                                 cylE5_r4 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,5));
tp@0:                                 cylE4_r5 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 6
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,5));
tp@0:                                 cylE6_r5 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,6));
tp@0:                                 cylE5_r6 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 7
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,6));
tp@0:                                 cylE7_r6 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,7));
tp@0:                                 cylE6_r7 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             if Ndifforder >= 8
tp@0:                                 index1 = reftoshortlistE(edgeshortlist(ii,8),edgeshortlist(ii,7));
tp@0:                                 cylE8_r7 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0:                                 index2 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,8));
tp@0:                                 cylE7_r8 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0:                             end
tp@0:                             
tp@0:                        else
tp@0:                             error('STRANGE TO END UP HERE??')
tp@0:                             if Ndifforder == 2
tp@0:                                 cylE1_r2 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))]; 
tp@0:                                 cylE2_r1 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))];
tp@0:                             else
tp@0:                                 error(['ERROR: Geometries with a single edge can not handle difforder >= 3'])
tp@0:                             end
tp@0:                        end
tp@0:                         
tp@0:                     for jj = 1:ncombs               
tp@0: 
tp@0:                         if multfactors(ivlonglist(iv(jj)))
tp@0:                             
tp@0:                             if Ndifforder == 2
tp@0:                                 
tp@0:                                 cylE1_r2frac = cylE1_r2;
tp@0:                                 e1length = edgelengthvec(edgeshortlist(ii,1));
tp@0:                                 cylE2_r1frac = cylE2_r1;
tp@0:                                 e2length = edgelengthvec(edgeshortlist(ii,2));
tp@0:                                 
tp@0:                                 if firstedgestart(jj) ~= 0
tp@0:                                     cylE1_r2frac(1,3) = cylE1_r2frac(1,3) + e1length*firstedgestart(jj);
tp@0:                                 end
tp@0:                                 if firstedgeend(jj) ~= 1
tp@0:                                     cylE1_r2frac(2,3) = cylE1_r2frac(1,3) + e1length*(1-firstedgeend(jj));
tp@0:                                 end
tp@0:                                 if lastedgestart(jj) ~= 0
tp@0:                                     cylE2_r1frac(1,3) = cylE2_r1frac(1,3) + e2length*lastedgestart(jj);
tp@0:                                 end
tp@0:                                 if lastedgeend(jj) ~= 1
tp@0:                                     cylE2_r1frac(2,3) = cylE2_r1frac(1,3)+ e2length*(1-lastedgeend(jj));
tp@0:                                 end
tp@0:              
tp@0:                                 if midspecpattern(iv(jj))~=0
tp@0:                                     
tp@0:                                     % For the cases with specular reflections
tp@0:                                     % in-between a double diffraction we must
tp@0:                                     % mirror the two edges to calculate the
tp@0:                                     % relative-to-edge cylindrical coordinates.
tp@0:                                     
tp@0:                                     nspec = sum(midspecpattern(iv(jj),:)>0);
tp@0: 			
tp@0:                                     % Mirror edge 2 in all the specular reflection
tp@0:                                     % planes, in reversed order
tp@0:                                     edgestartmirr = edgestartcoords(edgeshortlist(ii,2),:);
tp@0:                                     edgeendmirr =   edgeendcoords(edgeshortlist(ii,2),:);
tp@0:                                     edgevector = edgeendmirr - edgestartmirr;
tp@0:                                     if lastedgestart(jj) ~= 0
tp@0:                                         edgestartmirr = edgestartmirr + edgevector*lastedgestart(jj);
tp@0:                                     else
tp@0:                                         edgestartmirr = edgestartmirr + edgevector*1e-6;
tp@0:                                     end
tp@0:                                     if lastedgeend(jj) ~= 1
tp@0:                                         edgeendmirr =   edgestartmirr + edgevector*(1-lastedgeend(jj));
tp@0:                                     else
tp@0:                                         edgeendmirr =   edgestartmirr + edgevector*(1-1e-6);
tp@0:                                     end
tp@0:                                     edgerefcoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)];
tp@0:                                     edgerefnvec =   edgenvecs(edgeshortlist(ii,1),:);
tp@0: 			
tp@0:                                     % If we have a specular reflection in a plane
tp@0:                                     % which is perpendicular to the edge plane, we
tp@0:                                     % should nudge the mirrored edge out a bit so
tp@0:                                     % that there is no 0/(2*pi) mistake
tp@0:                                     if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) )
tp@0:                                         vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr;
tp@0:                                         edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10;
tp@0:                                         vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr;
tp@0:                                         edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10;
tp@0:                                     end
tp@0:                                     xis = [edgestartmirr;edgeendmirr];
tp@0:                                     for kk = nspec:-1:1
tp@0:                                         xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1);
tp@0:                                     end
tp@0:                                     [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec);
tp@0:                                     cylE2mirr_r1 = [rstart thetastart zstart;rend thetaend zend];
tp@0: 			
tp@0:                                     % Mirror edge 1 in all the specular reflection
tp@0:                                     % planes, in forward order
tp@0:                                     edgestartmirr = edgestartcoords(edgeshortlist(ii,1),:);
tp@0:                                     edgeendmirr =   edgeendcoords(edgeshortlist(ii,1),:);
tp@0:                                     edgevector = edgeendmirr - edgestartmirr;
tp@0:                                     if firstedgestart(jj) ~= 0
tp@0:                                         edgestartmirr = edgestartmirr + edgevector*firstedgestart(jj);
tp@0:                                     else
tp@0:                                         edgestartmirr = edgestartmirr + edgevector*1e-6;
tp@0:                                     end
tp@0:                                     if firstedgeend(jj) ~= 1
tp@0:                                         edgeendmirr =   edgestartmirr + edgevector*(1-firstedgeend(jj));
tp@0:                                     else
tp@0:                                         edgeendmirr =   edgestartmirr + edgevector*(1-1e-6);                               
tp@0:                                     end
tp@0:                                     edgerefcoords = [edgestartcoords(edgeshortlist(ii,2),:);edgeendcoords(edgeshortlist(ii,2),:)];
tp@0:                                     edgerefnvec =   edgenvecs(edgeshortlist(ii,2),:);
tp@0: 			
tp@0:                                     % Normally, when there is a specular reflection
tp@0:                                     % in-between, the diffraction path will not be
tp@0:                                     % along a plane, unless: 
tp@0:                                     % 1. The same edge is involved twice, with a
tp@0:                                     % reflection in a perpendicular plane
tp@0:                                     % in-between.
tp@0:                                     % 2. See below
tp@0:                                     pathalongplanewmidspec = 0;
tp@0:                                     if edgeshortlist(ii,1) == edgeshortlist(ii,2)
tp@0:                                         if thetastart == 0 | thetastart == (2*pi-closwedangvec(edgeshortlist(ii,2)))
tp@0:                                             pathalongplanewmidspec = 1;
tp@0:                                         end
tp@0:                                     end
tp@0:                                     
tp@0:                                     % If we have a specular reflection in a plane
tp@0:                                     % which is perpendicular to the edge plane, we
tp@0:                                     % should nudge the mirrored edge out a bit so
tp@0:                                     % that there is no 0/(2*pi) mistake
tp@0:                                     %
tp@0:                                     % We also have case 2 here (see above) for when
tp@0:                                     % we could have a pathalongplanewmidspec
tp@0:                                     % 2. When two different edges have a perpendicular reflection
tp@0:                                     % in-between 
tp@0: 			
tp@0:                                     if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) )
tp@0:                                         vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr;
tp@0:                                         edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10;
tp@0:                                         vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr;
tp@0:                                         edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10;
tp@0:                                         pathalongplanewmidspec = 1;
tp@0:                                     end
tp@0:                                     xis = [edgestartmirr;edgeendmirr];
tp@0:                                     for kk = 1:nspec
tp@0:                                         xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1);
tp@0:                                     end
tp@0:                                     [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec);
tp@0:                                     cylE1mirr_r2 = [rstart thetastart zstart;rend thetaend zend];
tp@0:                                     
tp@0:                                     [irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2mirr_r1,cylE1mirr_r2,...
tp@0:                                         nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),...
tp@0:                                         EDcalcmethod,pathalongplanewmidspec,Rstart,bc,CAIR,FSAMP);
tp@0:                                     
tp@0:                                 else   %  ....   if midspecpattern(iv(jj))~=0
tp@0: 			
tp@0: 									[irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2_r1frac,cylE1_r2frac,...
tp@0: 										nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),...
tp@0: 										EDcalcmethod,pathalongplane,Rstart,bc,CAIR,FSAMP);
tp@0:                                 end   %  ....   if midspecpattern(iv(jj))~=0
tp@0: 			
tp@0:                                 if SHOWTEXT >= 6
tp@0:                                     figure(1)
tp@0:                                     plot(irnew)
tp@0:                                    figure(2)
tp@0:                                     plot(irdiff)
tp@0:                                     pause
tp@0:                                 end  
tp@0:                                 IRDIFFVEC = [IRDIFFVEC;sum(irnew)];
tp@0:                                 			                
tp@0:                             elseif Ndifforder == 3,    %   if Ndifforder == 2			
tp@0: 			
tp@0: 								for kk = 1:newndivvec(1)
tp@0: 									if SHOWTEXT >= 5
tp@0: 										disp(['   ',int2str(kk),' of ',int2str(newndivvec(1))]) 
tp@0: 									end
tp@0: 									BIGEDGE1stvalue = (kk-0.5)./newndivvec(1);
tp@0: 									wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3];
tp@0: 									
tp@0: 									[irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,...
tp@0: 										nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',...
tp@0: 										dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue);
tp@0: 									irnewpartition = real(irnewpartition);
tp@0: 																			
tp@0: 									if kk == 1
tp@0: 										irnew = irnewpartition;    
tp@0: 									else
tp@0: 										lengthaddition = length(irnewpartition);
tp@0: 										lengthaccum = length(irnew);
tp@0: 										if lengthaddition > lengthaccum
tp@0: 											irnew = [irnew;zeros(lengthaddition-lengthaccum,1)];    
tp@0: 										end
tp@0: 										irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition; 
tp@0: 									end
tp@0: 								end
tp@0: 			
tp@0:                                 if SHOWTEXT >= 6
tp@0:                                     sum(irnew)
tp@0:                                     plot(irnew)
tp@0:                                     pause
tp@0:                                 end
tp@0:                             end
tp@0:                             if Ndifforder == 4
tp@0:                                     wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4];                                
tp@0:                             elseif Ndifforder == 5
tp@0:                                     wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5];
tp@0:                             elseif Ndifforder == 6
tp@0:                                     wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6];
tp@0:                             elseif Ndifforder == 7
tp@0:                                     wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7];
tp@0:                             elseif Ndifforder == 8,                                
tp@0:                                     wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7;cylE8_r7;cylE7_r8];
tp@0:                             end
tp@0:                             
tp@0:                             if Ndifforder >= 4,    
tp@0:                                 
tp@0:                                 for kk = 1:newndivvec(1)
tp@0:                                     if SHOWTEXT >= 5
tp@0:                                         disp(['   ',int2str(kk),' of ',int2str(newndivvec(1))]) 
tp@0:                                     end
tp@0:                                     BIGEDGE1stvalue = (kk-0.5)./newndivvec(1);
tp@0: 
tp@0:                                     [irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,...
tp@0:                                         nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',...
tp@0:                                         dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue);
tp@0:                                     irnewpartition = real(irnewpartition);
tp@0:                                         
tp@0:                                     if kk == 1
tp@0:                                         irnew = irnewpartition;
tp@0:                                     else
tp@0:                                         lengthaddition = length(irnewpartition);
tp@0:                                         lengthaccum = length(irnew);
tp@0:                                         if lengthaddition > lengthaccum
tp@0:                                             irnew = [irnew;zeros(lengthaddition-lengthaccum,1)];    
tp@0:                                         end
tp@0:                                         irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition; 
tp@0:                                     end
tp@0:                                 end
tp@0:                             end
tp@0:                                                             
tp@0:                             if SHOWTEXT >= 5
tp@0:                                 varname = ['allirs',int2str(Ndifforder)];
tp@0:                                 if exist(varname) == 0
tp@0:                                         eval([varname,' = irnew;'])
tp@0:                                 else
tp@0:                                         lengthaddition = length(irnew);
tp@0:                                         eval(['lengthaccum = size(',varname,',1);'])
tp@0:                                         if lengthaddition > lengthaccum
tp@0: 		                                    eval([varname,' = [',varname,';zeros(lengthaddition-lengthaccum,size(',varname,',2))];'])
tp@0:                                         end
tp@0:                                         eval([varname,' = [',varname,' zeros(size(',varname,',1),1)];'])
tp@0:                                         eval([varname,'(1:length(irnew),end) = irnew;'])
tp@0:                                 end
tp@0:                              end
tp@0:                             
tp@0:                             % Decide whether the IR should be boosted or not
tp@0: 			
tp@0:                             boostfactor = 1;
tp@0:                             if souspecboost ~= 1 | recspecboost ~= 1
tp@0:                             
tp@0:                                 if prespecpresent(iv(jj)) == 1
tp@0:                                     boostfactor = souspecboost;    
tp@0:                                 else
tp@0:                                     if ~isempty(Sinsideplanenumber)
tp@0:                                         if Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),1) & Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),2), 
tp@0:                                             boostfactor = souspecboost;
tp@0:                                         end
tp@0:                                     end    
tp@0: 			
tp@0:                                 end
tp@0:                                 if postspecpresent(iv(jj)) == 1
tp@0:                                     boostfactor = boostfactor*recspecboost;    
tp@0:                                 else
tp@0:                                     if ~isempty(Rinsideplanenumber)
tp@0:                                         if Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),1) & Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),2), 
tp@0:                                             boostfactor = boostfactor*recspecboost;
tp@0:                                         end
tp@0:                                     end
tp@0:                                 end
tp@0:                             
tp@0:                             end   % ...   if souspecboost ~= 1 | recspecboost ~= 1
tp@0:                             
tp@0:                             % For thin plates, we must have a boost factor!
tp@0:                             % This is because there will be multiple equivalent
tp@0:                             % combinations passing on the rear side of the thin plate
tp@0:                             
tp@0:                             if all( nyvec(edgeshortlist(ii,:)) == 0.5 )
tp@0:                                 boostfactor = boostfactor*2^(Ndifforder-1);                        
tp@0:                             end
tp@0: 			
tp@0:                             if boostfactor ~= 1
tp@0:                                 irnew = irnew*boostfactor;    
tp@0:                             end
tp@0:                             
tp@0:                             irnew = irnew*double(multfactors(ivlonglist(iv(jj))));
tp@0: 			
tp@0:                             [ndiff,ncolsdiff] = size(irdiff);
tp@0:                             nnew = size(irnew,1);
tp@0:                             if saveindividualdiffirs(1) == 0
tp@0:                                 if nnew > ndiff
tp@0:                                    irdiff = [irdiff;zeros(nnew-ndiff,1)];
tp@0:                                 end
tp@0:                                 irdiff(1:nnew) = irdiff(1:nnew) + irnew;
tp@0:                             else
tp@0:                                 if Ndifforder > ncolsdiff
tp@0:                                     irdiff = [irdiff zeros(ndiff,Ndifforder-ncolsdiff)];    
tp@0:                                     ncolsdiff = ncolsdiff+1;
tp@0:                                 end
tp@0:                                 if nnew > ndiff
tp@0:                                    irdiff = [irdiff;zeros(nnew-ndiff,Ndifforder)];
tp@0:                                 end
tp@0:                                 irdiff(1:nnew,Ndifforder) = irdiff(1:nnew,Ndifforder) + irnew;
tp@0:                             end
tp@0:                             
tp@0:                         end    % ...  if multfactors(ivlonglist(iv(jj)))
tp@0:                         
tp@0:                     end    % ....for jj = 1:ncombs
tp@0: 
tp@0:                 else   %  ...    if any(multfactors(ivlonglist(iv)))
tp@0:                     if SHOWTEXT >= 4
tp@0:                         disp(['      Combination not computed because of repetition'])                         
tp@0:                     end
tp@0:                     
tp@0:                 end    % ...    if any(multfactors(ivlonglist(iv)))
tp@0:                 
tp@0:             end   % .... for ii = 1: length(ivec)
tp@0:         
tp@0:         end   % ...     if any(difforderinlist==Ndifforder)
tp@0:     
tp@0:         if size(irdiff,2) < Ndifforder & saveindividualdiffirs(1) == 1
tp@0:             irdiff = [irdiff zeros(size(irdiff,1),Ndifforder-size(irdiff,2))];    
tp@0:         end
tp@0: 	end    % ... for Ndifforder = 2:specorder
tp@0: end     % .... if userwantsdiff2 == 1
tp@0: 
tp@0: ntot = length(irtot);
tp@0: ndirect = length(irdirect);
tp@0: ngeom = size(irgeom,1);
tp@0: ndiff = size(irdiff,1);
tp@0: 
tp@0: if ndirect > ntot
tp@0:    irtot = [irtot;zeros(ndirect-ntot,1)];
tp@0:     ntot = length(irtot);
tp@0: end
tp@0: irtot(1:ndirect) = irtot(1:ndirect) + irdirect;
tp@0: 
tp@0: if ngeom > ntot
tp@0:    irtot = [irtot;zeros(ngeom-ntot,1)];
tp@0:     ntot = length(irtot);
tp@0: end
tp@0: irtot(1:ngeom) = irtot(1:ngeom) + irgeom;
tp@0: 
tp@0: if ndiff > ntot
tp@0:    irtot = [irtot;zeros(ndiff-ntot,1)];
tp@0:     ntot = length(irtot);
tp@0: end
tp@0: if saveindividualdiffirs(1) == 0 | userwantsdiff2 == 0
tp@0:     irtot(1:ndiff) = irtot(1:ndiff) + irdiff;
tp@0: else
tp@0:      irtot(1:ndiff) = irtot(1:ndiff) + sum(irdiff.').';    
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Make the IRs the same length
tp@0: 
tp@0: nmax = max([ length(irtot) size(irgeom,1) size(irdiff,1) length(irdirect)]);
tp@0: 
tp@0: if length(irtot) < nmax
tp@0: 	irtot = [irtot;zeros(nmax-length(irtot),1)];
tp@0: end
tp@0: if length(irdirect) < nmax
tp@0: 	irdirect = [irdirect;zeros(nmax-length(irdirect),1)];
tp@0: end
tp@0: if length(irgeom) < nmax
tp@0: 	irgeom = [irgeom;zeros(nmax-size(irgeom,1),size(irgeom,2))];
tp@0: end
tp@0: if length(irdiff) < nmax
tp@0: 	irdiff = [irdiff;zeros(nmax-size(irdiff,1),size(irdiff,2))];
tp@0: end
tp@0: 
tp@0: %-------------------------------------------------------
tp@0: % Save the IRs
tp@0: 
tp@0: irtot = sparse(irtot);
tp@0: irgeom = sparse(irgeom);
tp@0: irdirect = sparse(irdirect);
tp@0: irdiff = sparse(irdiff);
tp@0: Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod'];
tp@0: 
tp@0: if length(saveindividualdiffirs) > 1
tp@0: 	if saveindividualdiffirs(2) == 1
tp@0: 		Varlist = [Varlist,' alldiffirs alldiffirsextradata'];
tp@0: 	end
tp@0: end
tp@0: 
tp@0: if SHOWTEXT >= 5
tp@0:     if specorder >= 2
tp@0:         if exist('allirs2') == 0
tp@0:             allirs2 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs2'];
tp@0:     end       
tp@0:     if specorder >= 3
tp@0:         if exist('allirs3') == 0
tp@0:             allirs3 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs3'];
tp@0:     end       
tp@0:     if specorder >= 4
tp@0:         if exist('allirs4') == 0
tp@0:             allirs4 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs4'];
tp@0:     end       
tp@0:     if specorder >= 5
tp@0:         if exist('allirs5') == 0
tp@0:             allirs5 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs5'];
tp@0:     end       
tp@0:     if specorder >= 6
tp@0:         if exist('allirs6') == 0
tp@0:             allirs6 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs6'];
tp@0:     end       
tp@0:     if specorder >= 7
tp@0:         if exist('allirs7') == 0
tp@0:             allirs7 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs7'];
tp@0:     end       
tp@0:     if specorder >= 8
tp@0:         if exist('allirs8') == 0
tp@0:             allirs8 = [];    
tp@0:         end
tp@0:         Varlist = [Varlist,' allirs8'];
tp@0:     end       
tp@0: end
tp@0: 
tp@0: eval(['save ',edirfile,Varlist])