tp@0: function EDB1main(EDsetupfile)
tp@0: % EDB1main - Calculates the specular and edge diffraction IRs and saves them in a file.
tp@0: % Calculates the specular and edge diffraction IRs and saves them in a file or
tp@0: % a series of files, if there are many sources and receivers.The
tp@0: % calculation goes through three stages:
tp@0: %   1. Geometrical precalculations: define edges and determine visibility
tp@0: %      between edges, planes, sources and receivers.
tp@0: %   2. Find the valid sound paths and save a list of them
tp@0: %   3. Go through the list and construct IRs
tp@0: %
tp@0: % This version, EDB1main, includes the option of dividing the
tp@0: % calculations into batches. When there are many receivers, the geometrical
tp@0: % precalculations might overload the available memory since all receivers
tp@0: % are treated simulataneously for each source.
tp@0: %
tp@0: % Input parameters:
tp@0: %	EDsetupfile 	(optional)	This file contains all parameter settings for the calculations.
tp@0: %					If no filename is given, a file open window will appear and the
tp@0: %					desired EDsetupfile can be chosen.
tp@0: %                   The file could either be a .m file (the "EDsetupfile")
tp@0: %                   or a .mat file (the "setupmatfile").
tp@0: % Input parameters in the EDsetupfile:
tp@0: %	FSAMP, CAIR, RHOAIR, SHOWTEXT
tp@0: %						(global)	The sampling frequency, speed of sound, density of air
tp@0: %									and the control of how much text should
tp@0: %									be plotted (SHOWTEXT could be 0,1,2,3,4).
tp@0: %   SUPPRESSFILES (global, optional) If this variable is introduced in the
tp@0: %                       setupfile as a global variable, with the value 1, then the files
tp@0: %                       'xxx_ISEStree.mat' will not be generated. Instead values are passed
tp@0: %                       directly.
tp@0: %	Filepath			The path to the folder where all output data files will be stored.
tp@0: %	Filestem			The first part of the output data files will be given this name
tp@0: %	CADfile			    The name of the .CAD file where the geometry is given, including path
tp@0: %   open_or_closed_model      Specifies whether the geometrical model is open or closed:
tp@0: %                       'o'     Open model
tp@0: %                       'c'     Closed model
tp@0: %   int_or_ext_model    Specifies whether the geometrical model is interior or exterior:
tp@0: %                       'i'     Interior model (NB! An interior model must be closed, 
tp@0: %                               otherwise it should be called an exterior model)
tp@0: %                       'e'     Exterior model
tp@0: %	EDcalcmethod	    Specifies the edge diffraction method. Must have one of the following values:
tp@0: %						'v'		Vanderkooy
tp@0: %						'n'		Svensson et al
tp@0: %						'k'		Kirchhoff approximation NB! The Kirchhoff
tp@0: %						        approximation gives only first-order diffraction.
tp@0: %	directsound		0 or 1
tp@0: %	specorder		The maximum order of specular reflections. (To get all possible combinations
tp@0: %                   specorder and difforder should be given the same value).
tp@0: %	difforder		The maximum order of diffraction components. (difforder can never be higher
tp@0: %                   than specorder, i.e., if you want difforder = 2, you must set specorder to 2 or higher)
tp@0: %   elemsize		A vector of size [1,difforder] which specifies the size of the edge subdivision
tp@0: %                   for edge diffraction calculations. A value of 0.2 means that the edge elements
tp@0: %                   are 5 times (1/0.2 =5) the wavelength at FSAMP.
tp@0: %				    A higher value gives more accurate edge diffraction but is also more time consuming. 
tp@0: %                   Larger and larger elements could be employed for higher-order diffraction.
tp@0: %                   Suggested values for high accuracy: elemsize = [1 0.5 0.25 0.125 ....].
tp@0: %   nedgesubs       All edges are subdivided very coarsely into edge segments, for visibility checks.
tp@0: %                   The value of nedgesubs states how many such segments each edge will be divided into.
tp@0: %                   The minimum value is 2. NB! All edges will be divided into equally many segments
tp@0: %                   in order to employ efficient Matlab vector-based processing.
tp@0: %	sources			A vector of the source coordinates, size [nsources,3], or source names/numbers
tp@0: %                   size [nsources,1] or [nsources,2] (if they should be taken from the CADfile.)
tp@0: %	receivers		A vector of the receiver coordinates or receiver numbers (if they should be taken
tp@0: %					from the CADfile. Same syntax as for sources.
tp@0: %   nSbatches,nRbatches
tp@0: %                   (optional) An integer >= 1 which defines how many chunks the number of sources/receivers
tp@0: %                   will be divided into. Default = 1. nSbatches/nRbatches = 1 means that all sources/receivers are 
tp@0: %                   treated simultaneously. nSbatches/nRbatches = 2 means that the list of sources/receivers
tp@0: %                   is divided into two chunks etc. It is fastest with nSbatches/nRbatches = 1 but it also
tp@0: %                   requires the largest amount of memory. (For large problems, there might even be no time
tp@0: %                   penalty in choosing nRbatches = the number of receivers!)
tp@0: %   soulist_for_findpaths/reclist_for_findpaths
tp@0: %                   (optional) A list of integers that specifies which of the sources/receivers that
tp@0: %                   should be run for the "find paths" stage. Default = [1:nsources]/[1:nreceivers].
tp@0: %                   This is useful for huge problems where one might want to stop (with CTRL-C) and restart calculations.
tp@0: %   soulist_for_calcirs/reclist_for_calcirs
tp@0: %                   (optional) A list of integers that specifies which of the sources/receivers that
tp@0: %                   should be run for the "calc irs" stage. Default = [1:nsources]/[1:nreceivers].
tp@0: %                   This is useful for huge problems where one might want to stop (with CTRL-C) and restart calculations.
tp@0: %	firstskipcorner (optional) All edges including at least one corner with this number or
tp@0: %             		higher will be excluded from the calculations. Default: 1e6
tp@0: %	Rstart 			(optional) The distance that should correspond to the time zero for
tp@0: %             		all IRs. Default = 0. NB! With Rstart = 0, all IRs will include the initial time
tp@0: %                   delay that is caused by the propagation time from the source to the receiver. For
tp@0: %             		long distances and high FSAMP, this can lead to many initial samples that are zero.
tp@0: %                   If you know that all source-receiver combinations have a certain minimum
tp@0: %             		distance, this minimum distance (or a value slightly below) can be given as Rstart.
tp@0: %                   If Rstart is given a value which is smaller than the shortest distance encountered, 
tp@0: %                   an error will occur.
tp@0: %
tp@0: % Optional input parameters in the EDsetupfile:
tp@0: %	cadgeofile	    If this parameter is specified, this file will be read rather than the CADfile.
tp@0: %                   Saves some time for big geometries.
tp@0: %	eddatafile	    If this parameter is specified, the function EDB1edgeo is not called. Saves
tp@0: %					some time.
tp@0: %	srdatafile	    If this parameter is specified, the function EDB1srgeo is not called. Saves
tp@0: %					some time.
tp@0: %   restartnumber   If this parameter is specified, the calculations will start at receiver number
tp@0: %                   restartnumber, rather than number 1. This is useful if a calculation for many
tp@0: %                   receivers can not be completed in one run.
tp@0: %   stopnumber      If this parameter is specified, the calculations will stop for receiver number
tp@0: %                   stopnumber.
tp@0: %   saveindividualdiffirs   A vector with two values, 0 or 1.
tp@0: %					If the first value is given the value 1, the diffraction
tp@0: %                   IR will have one column for each order. Default is that
tp@0: %                   all orders are summed into a single irdiff.
tp@0: %					If the second value is given the value 1, all individual diffraction irs
tp@0: %					will be saved on individual lines in a large matrix called 'alldiffirs'.
tp@0: %				 	This matrix has only single precision.
tp@0: %					Default value: [0 0]
tp@0: %   doaddsources    0 or 1, (default value 0) indicating whether:
tp@0: %                       0: every individual S-to-R IR should be saved
tp@0: %                       1: a single IR is saved for each R, which is
tp@0: %                       constructed by first computing all individual
tp@0: %                       S-to-R IRs, saving them in a temp file, and then
tp@0: %                       add IRs from all sources together. This way a
tp@0: %                       loudspeaker element can be simulated by a
tp@0: %                       distribution of point sources.
tp@0: %   dobackscatter   0 or 1, (default value 0), indicating whether:
tp@0: %                       0: every individual S-to-R IR should be saved
tp@0: %                       1: only S1 to R1, S2 to R2, S3 to R3 combos
tp@0: %                       are calculated.
tp@0: %
tp@0: % Output parameters:
tp@0: % 	The following IRs are calculated and saved in the output file(s). Any of them can be
tp@0: % 	excluded from the calculations by setting the parameters in the setup file.
tp@0: %
tp@0: %	irtot    		the total impulse response
tp@0: %	irdirect 		the direct sound
tp@0: %  	irgeom   		the purely geometrical acoustics components, i.e. specular reflections.
tp@0: %  	irdiff   		the diffracted components. This will innclude all combinations of specular and  
tp@0: %                   diffraction, and first-order as well as higher-order combinations.
tp@0: %
tp@0: % Note that these variables, together with those in the EDsetupfile are saved in a results file
tp@0: %	which gets the name in Filestem (see below) + '_SS_RR_ir.mat' (if addsources = 0) or
tp@0: %	'_RR_ir.mat' (if addsources = 1).
tp@0: %
tp@0: % Note also that the impulse responses are scaled such that the direct sound gets an amplitude of
tp@0: % 1/distance. This means that the impulse response represents a system with the sound pressure as
tp@0: % output signal and an input signal = RHOAIR*Volumeacc(t)/(4*pi) where Volumeacc(t) is
tp@0: % the volume acceleration of the monopole source.
tp@0: %
tp@0: % A logfile containing the computation time for the three stages is
tp@0: % created, with the name Filestem + '_log.mat'
tp@0: %
tp@0: % Uses functions EDB1readcad EDB1readsetup EDB1edgeo
tp@0: % EDB1Sgeo EDB1Rgeo EDB1findISEStree EDB1findpathsISES EDB1makeirs
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) 20100211
tp@0: %
tp@0: % EDB1main(EDsetupfile);
tp@0: 
tp@0: global FSAMP CAIR RHOAIR SHOWTEXT JJ JJnumbofchars SUPPRESSFILES
tp@0: global POTENTIALISES ISCOORDS IVNDIFFMATRIX
tp@0: global IVNSPECMATRIX ORIGINSFROM ISESVISIBILITY REFLORDER
tp@0: 
tp@0: global IRDIFFVEC
tp@0: 
tp@0: FSAMP = []; CAIR = []; RHOAIR = [];  SHOWTEXT = [];
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % First run the setup m file, just to get access to the CAD-file name
tp@0: % so that the CAD-file can be converted into a cadgeofile.
tp@0: % (This conversion must be done before the setup m-file is converted
tp@0: % to a setup mat-file since the source and receiver coordinates should
tp@0: % be specified in the setupmatfile, and the sources and receivers might
tp@0: % have to be read from the CAD-file).
tp@0: %
tp@0: % First we make a little check of a common mistake: check if the paths
tp@0: % seem to be for the correct computer type.
tp@0: % The EDsetupfile-path does not need to be checked if a file window is
tp@0: % presented. The CAD-file path always needs to be checked.
tp@0: 
tp@0: % We check if the type of setupfile:
tp@0: %   .m    setupfiletype = 1;
tp@0: %   .mat  setupfiletype = 2;
tp@0: 
tp@0: setupfiletype = 1;
tp@0: 
tp@0: if nargin == 0
tp@0: 	[EDsetupfile,EDsetupfilepath] = uigetfile('*.m','Please select the EDsetupfile');
tp@0:     [temp1,EDsetupfile,temp2] = fileparts([EDsetupfilepath,EDsetupfile]);
tp@0: 
tp@0:     EDsetupfile = [EDsetupfilepath,EDsetupfile];
tp@0:     
tp@0: 	if ~isstr(EDsetupfile) | isempty(EDsetupfile)
tp@0: 		return
tp@0: 	end
tp@0: else
tp@0:     [EDsetupfilepath,tempfilestem,fileext] = fileparts(EDsetupfile);
tp@0:     EDsetupfilepath = [EDsetupfilepath,filesep];
tp@0:     if length(fileext) == 4
tp@0:         if fileext == '.mat'
tp@0:             setupfiletype = 2; 
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
luis@9: % need to read setup file manually --luisf.
tp@0: 
luis@9: % if setupfiletype == 1
luis@9: %     comptype = computer;
luis@9: %     if length(comptype) == 3
luis@9: %         comptype = [comptype,'P'];
luis@9: %     end
luis@9: %     if comptype(1:4) == 'MACI', 
luis@9: %         [temp1,tempfilestem,temp2] = fileparts(EDsetupfile);        
luis@9: %         eval([tempfilestem])                       
luis@9: %     else    
luis@9: %         eval(['run ''',EDsetupfile,''''])
luis@9: %     end
luis@9: % else
luis@9: %     eval(['load ''',EDsetupfile,''''])
luis@9: % end
luis@9: 
luis@9: eval(char(textread(EDsetupfile, '%s', 'whitespace', '')));
luis@9: % ~luisf
tp@0: 
tp@0: if SHOWTEXT >= 1
tp@0: 	disp(' ')
tp@0: 	disp('####################################################################')
tp@0: 	disp('#')
tp@0: 	disp('#  EDB1mainISESx, version 11 February 2010')
tp@0: 	disp('#')
tp@0: end
tp@0: 
tp@0: if SHOWTEXT >= 2
tp@0: 	disp('#')
tp@0: 	disp('# This version does the following:')
tp@0: 	disp('# *  Caclulates impulse responses that include specular reflections up to any order and') 
tp@0:     disp('#    multiple diffraction up to order 6. Combinations are also taken into account, such as')
tp@0:     disp('#    specular-diffraction-specular but higher-order diffractions must be')
tp@0:     disp('#    in a single sequence, possibly with specular reflections before and after.')
tp@0:     disp('#    There will be memory problems for higher-order specular reflections due to a')
tp@0:     disp('#    memory-hungry vectorized version of the image source method.')
tp@0:     disp('# *  Combinations with two diffractions that have specular reflections inbetween')
tp@0:     disp('#    are handled, but not with any specular reflections after the last diffraction.')
tp@0: 	disp('# *  Partly visible edges are handled in first-order diffraction')
tp@0: 	disp('#    calculations. The accuracy depends on how high value of')
tp@0: 	disp('#    nedgesubs that is chosen. Edges that are split into two or more subedges are handled properly.')
tp@0:     disp('# *  Specular reflections in the IR are handled by splitting a pulse between two consecutive time slots.')
tp@0:     disp('#    This gives quite a severe low-pass filter roll-off, but zero phase-error.')
tp@0:     disp('#    Choose a high over-sampling ratio to minimize this magnitude error.')
tp@0:     disp('#    In addition, the low-pass filter effect can be avoided altogether for the direct sound by setting')
tp@0:     disp('#    Rstart to exactly the source-receiver distance. Thereby the direct sound pulse will occur exactly in the')
tp@0:     disp('#    first time sample. For more tweaking, the three values Rstart and CAIR and FSAMP should be possible')
tp@0:     disp('#    to adjust so that one, two or three geometrical pulses arrive exactly at sample times.')
tp@0:     disp('# *  An accurate integration technique is used for first-order diffraction so that receiver positions close')
tp@0:     disp('#    to zone boundaries are handled without problem. Exactly at the zone boundary the specular reflection')
tp@0:     disp('#    or direct sound gets half the amplitude (See below though).')
tp@0: 	disp('#')
tp@0: 	disp('# Limitations with this version:')
tp@0: 	disp('# 1. Some, but not all, possible combinations where two or more diffracted')
tp@0: 	disp('#    reflections have specular reflections in-between are calculated. ')
tp@0:     disp('#    This requires a bit of work to complete.')
tp@0: 	disp('# 3. In-plane visibility of edges is only partly checked, i.e., if a plane has indents')
tp@0: 	disp('#    so that an edge can not see all other edges of its own plane.')
tp@0: 	disp('# 4. Part-visibility of higher-order diffraction is not checked. These components are treated')
tp@0:     disp('#    as fully visible or totally invisible. However, if the first or last edge in the')
tp@0:     disp('#    sequence is partly visible, this will be used.')
tp@0: 	disp('# 5. Non-rigid specular reflections or edge diffraction for non-rigid planes are not implemented. ')
tp@0: 	disp('#    Pressure-release surfaces would be quite easy to implement.')
tp@0:     disp('#    Other surface impedances could be implemented since the hitpoints for all specular')
tp@0:     disp('#    reflections are stored in the edpaths output file.')
tp@0:     disp('# 6. Other receivers than omnidirectional microphones are not implemented. HRTFs or directional microphones')
tp@0:     disp('#    could be implemented since the hitpoints for all specular reflections are stored in the edpaths')
tp@0:     disp('#    so the incidence angles are easy to calculate. Reflections that end with a diffraction need some')
tp@0:     disp('#    more consideration.')
tp@0:     disp('# 7. The special cases where the receiver is exactly at a zone boundary is handled correctly only')
tp@0:     disp('#    for first-order diffraction. This takes a bit of work to implement for higher-order diffraction.')
tp@0:     disp('# 8. Sources and receivers can not be place directly at a surface. They must be raised by 1 mm or so.')
tp@0: 	disp('#')
tp@0: 	disp('# This version could be made more efficient by:')
tp@0: 	disp('# 1. Employing ray-tracing or beam-tracing for finding the visibility of')
tp@0: 	disp('#    plane to plane etc. Now, a simple plane-to-plane visibility check')
tp@0: 	disp('#    is performed: are planes in front of each other? The final visibility')
tp@0: 	disp('#    check is done for each image source to each receiver. An improved method would be more efficient')
tp@0: 	disp('#    for geometries with lots of indents and obscuring planes, e.g., as in a city')
tp@0: 	disp('#    geometry.')
tp@0: 	disp('# 2. Developing an algorithm that calculates the diffraction IR in an iterative process')
tp@0: 	disp('#    instead of restarting for each new order of diffraction.')
tp@0:     disp('# 3. The direct sound visibility test could be done earlier (in EDB1SRgeo).')
tp@0:     disp('# 4. Calculation of Image Receiver coordinates could probably be made more efficiently.')
tp@0:     disp('#')
tp@0: end
tp@0: 
tp@0: if SHOWTEXT >= 1
tp@0:     disp(' ')
tp@0:     disp(['   Using the setupfile: ',EDsetupfile])
tp@0:     disp(' ')
tp@0: end
tp@0: 
tp@0: % #########################################################################
tp@0: %   
tp@0: %      Geometry pre-calculations
tp@0: %   
tp@0: % #########################################################################
tp@0: 
tp@0: %---------------------------------------------------------------------------------------------
tp@0: % Convert the CAD-file into a mat file - or read the cadgeofile if it is specified
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: if SHOWTEXT >= 1
tp@0: 	disp(' ')
tp@0: 	disp('####################################################################')
tp@0: 	disp('#')
tp@0: 	disp('#  Geometry pre-calculations')
tp@0: 	disp(' ')
tp@0: end
tp@0: 
tp@0: if exist('cadgeofile') ~= 1
tp@0:     desiredname = [Filepath,Filestem,'_cadgeo'];
tp@0:     if SHOWTEXT >= 2
tp@0: 		disp(['   Reading CAD-file: ',CADfile,' and converting it to a cadgeofile: ',desiredname])
tp@0:     end
tp@0:     % Edited by APO: Added support for .ac files (only geometry)
tp@0:     [temp1,temp2,fileext] = fileparts(CADfile);
tp@0:     if(strcmp(fileext,'.ac'))
tp@0:         cadgeofile = EDB1readac(CADfile,desiredname);
tp@0:     else
tp@0:         cadgeofile = EDB1readcad(CADfile,desiredname);
tp@0:     end
tp@0: else
tp@0: 	if SHOWTEXT >= 2
tp@0: 		disp(['   Using existing cadgeofile: ',cadgeofile])
tp@0: 	end
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % Convert the EDsetupfile into a setupmatfile. If sources and receivers should
tp@0: % be read in the CADfile, add these to the setupmatfile.
tp@0: %
tp@0: % This is needed only if the setupfile was of type 1, that is, a .m file.
tp@0: 
tp@0: if setupfiletype == 1
tp@0:     setupmatfile = EDB1readsetup(EDsetupfile,cadgeofile);
tp@0: 
tp@0:     eval(['load ',setupmatfile])
tp@0: end
tp@0: 
tp@0: tcadgeo = etime(clock,t00);
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % Derive the edge parameters.
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: if calcpaths == 1
tp@0:     
tp@0:     if exist('eddatafile') ~= 1
tp@0: 	    desiredfile = [Filepath,Filestem,'_eddata'];
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Calculating edges and creating an eddatafile: ',desiredfile])
tp@0: 		end
tp@0: 		eddatafile = EDB1edgeox(cadgeofile,desiredfile,lower(open_or_closed_model(1)),lower(int_or_ext_model(1)),...
tp@0:             specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy);
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing eddatafile: ',eddatafile])
tp@0: 		end
tp@0: 	end
tp@0: else
tp@0: 	if exist('eddatafile') ~= 1
tp@0: 	    desiredfile = [Filepath,Filestem,'_eddata'];
tp@0:     	eddatafile = desiredfile;
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing eddatafile: ',eddatafile])
tp@0: 		end
tp@0: 	end
tp@0: end
tp@0: 
tp@0: tedgeo = etime(clock,t00);
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % If batches should be used, divide the list of receivers into batches.
tp@0: 
tp@0: nsources = size(sources,1);
tp@0: nreceivers = size(receivers,1);
tp@0: 
tp@0: if nSbatches == 1
tp@0:     doSbatches = 0;    
tp@0: elseif nSbatches > 1
tp@0:     if nSbatches > nsources
tp@0:         error('ERROR: nSbatches can not be larger than the number of sources')    
tp@0:     end
tp@0:     doSbatches = 1;
tp@0:     Sbatchlist = zeros(nSbatches,2);
tp@0:     nperbatch = ceil(nsources/nSbatches);
tp@0:     Sbatchlist(:,2) = [1:nSbatches].'*nperbatch;
tp@0:     Sbatchlist(:,1) = [0:nSbatches-1].'*nperbatch+1;
tp@0:     Sbatchlist(nSbatches,2) = nsources;
tp@0:     reftoSbatchlist = [1:nsources].';
tp@0:     reftoSbatchlist = reftoSbatchlist(:,ones(1,nSbatches));
tp@0:     compvec = Sbatchlist(:,1).';
tp@0:     reftoSbatchlist = reftoSbatchlist >= compvec(ones(nsources,1),:);
tp@0:     reftoSbatchlist = sum(reftoSbatchlist.').';
tp@0: end
tp@0: 
tp@0: if nRbatches == 1
tp@0:     doRbatches = 0;    
tp@0: elseif nRbatches > 1
tp@0:     if nRbatches > nreceivers
tp@0:         error('ERROR: nRbatches can not be larger than the number of receivers')    
tp@0:     end
tp@0:     doRbatches = 1;
tp@0:     Rbatchlist = zeros(nRbatches,2);
tp@0:     nperbatch = ceil(nreceivers/nRbatches);
tp@0:     Rbatchlist(:,2) = [1:nRbatches].'*nperbatch;
tp@0:     Rbatchlist(:,1) = [0:nRbatches-1].'*nperbatch+1;
tp@0:     Rbatchlist(nRbatches,2) = nreceivers;
tp@0:     reftoRbatchlist = [1:nreceivers].';
tp@0:     reftoRbatchlist = reftoRbatchlist(:,ones(1,nRbatches));
tp@0:     compvec = Rbatchlist(:,1).';
tp@0:     reftoRbatchlist = reftoRbatchlist >= compvec(ones(nreceivers,1),:);
tp@0:     reftoRbatchlist = sum(reftoRbatchlist.').';
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------------------------------------
tp@0: % Find out the S and R related parameters, make obstruction checks etc by calling 
tp@0: % EDB1Sgeo/EDB1Rgeo or load existing sdata/rdata files.
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: if calcpaths == 1
tp@0: 	if exist('sdatafile') ~= 1
tp@0:         if doSbatches == 0
tp@0: 		    desiredname = [Filepath,Filestem,'_sdata'];
tp@0: 			if SHOWTEXT >= 2
tp@0: 				disp(['   Calculating S parameters and creating an sdatafile: ',desiredname])
tp@0: 			end
tp@0: 
tp@0:             sdatafile = EDB1SorRgeo(eddatafile,desiredname,sources,'S',difforder,nedgesubs);
tp@0:         else
tp@0:             error(['ERROR: S-batches not implemented yet'])
tp@0:             for ii = 1:nSbatches
tp@0: 			    desiredname = [Filepath,Filestem,'_sdata'];
tp@0: 				if SHOWTEXT >= 2
tp@0: 					disp(['   Calculating S parameters and creating an sdatafile: ',desiredname,'_',int2str(ii)])
tp@0: 				end
tp@0:         		sdatafile = EDB1SorRgeo(eddatafile,[desiredname,'_',int2str(ii)],sources(Sbatchlist(ii,1):Sbatchlist(ii,2),:),'S',difforder,nedgesubs);
tp@0:                 sdatafile = desiredname;   % To trim off the _1 or _2 or...    
tp@0:             end
tp@0:         end
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing sdatafile: ',sdatafile])
tp@0: 		end
tp@0: 	end
tp@0: else
tp@0: 	if exist('sdatafile') ~= 1
tp@0: 	    desiredname = [Filepath,Filestem,'_sdata'];
tp@0:     	sdatafile = desiredname;
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing sdatafile: ',sdatafile])
tp@0: 		end
tp@0: 	end
tp@0: end
tp@0: 
tp@0: tsgeo = etime(clock,t00);
tp@0: 
tp@0: %------------------------------------
tp@0: % We check that there is not a small geometrical mistake which makes the
tp@0: % source be behind all planes - but only if the user has set SHOWTEXT.
tp@0: 
tp@0: if SHOWTEXT >= 2
tp@0:     eval(['load ',sdatafile])
tp@0:     if size(visplanesfroms,2) == 1
tp@0:         if ~any(visplanesfroms)
tp@0:             disp(' ')
tp@0:             disp('WARNING!!! The source(s) can not see a single plane!')
tp@0:             disp('           Check if this is correct! ')
tp@0:             disp('           A common cause is that the source is placed very close to a plane,')
tp@0:             disp('           but behind it. (CR => continue calculations)')
tp@0:             disp(' ')
tp@0:             pause
tp@0:         end                
tp@0:     else
tp@0:         iv = find(sum(visplanesfroms) == 0);
tp@0:         if ~isempty(iv)
tp@0:             disp(' ')
tp@0:             disp('WARNING!!! Some of the source(s) can not see a single plane!')
tp@0:             disp('           These are the sources number:')
tp@0:            disp(['           ',int2str(iv(:).')])
tp@0:             disp('           Check if this is correct! ')
tp@0:             disp('           A common cause is that the source is placed very close to a plane,')
tp@0:             disp('           but behind it. (CR => continue calculations)')
tp@0:             disp(' ')
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: if calcpaths == 1
tp@0: 	if exist('rdatafile') ~= 1
tp@0:         if doRbatches == 0
tp@0: 		    desiredname = [Filepath,Filestem,'_rdata'];
tp@0: 			if SHOWTEXT >= 2
tp@0: 				disp(['   Calculating R parameters and creating an rdatafile: ',desiredname])
tp@0: 			end
tp@0:     		rdatafile = EDB1SorRgeo(eddatafile,desiredname,receivers,'R',difforder,nedgesubs);
tp@0:         else
tp@0:             for ii = 1:nRbatches
tp@0: 			    desiredname = [Filepath,Filestem,'_rdata'];
tp@0: 				if SHOWTEXT >= 2
tp@0: 					disp(['   Calculating R parameters and creating an rdatafile: ',desiredname,'_',int2str(ii)])
tp@0: 				end
tp@0:                 
tp@0:                 rdatafile = EDB1SorRgeo(eddatafile,[desiredname,'_',int2str(ii)],receivers(Rbatchlist(ii,1):Rbatchlist(ii,2),:),'R',difforder,nedgesubs);
tp@0:                 rdatafile = desiredname;   % To trim off the _1 or _2 or...    
tp@0:             end
tp@0:         end
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing rdatafile: ',rdatafile])
tp@0: 		end
tp@0: 	end
tp@0: else
tp@0: 	if exist('rdatafile') ~= 1
tp@0: 	    desiredname = [Filepath,Filestem,'_rdata'];
tp@0:     	rdatafile = desiredname;
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing rdatafile: ',rdatafile])
tp@0: 		end
tp@0: 	end
tp@0: end
tp@0: 
tp@0: trgeo = etime(clock,t00);
tp@0: 
tp@0: %------------------------------------
tp@0: % We check that there is not a small geometrical mistake which makes the
tp@0: % receiver(s) be behind all planes - but only if the user has set SHOWTEXT.
tp@0: 
tp@0: if SHOWTEXT >= 2
tp@0:     if exist('nRbatches') == 0
tp@0:         nRbatches = 1;
tp@0:     end
tp@0:     if nRbatches == 1
tp@0:         eval(['load ',rdatafile])
tp@0:     else
tp@0:         [temp1,rdatafilestem,temp2] = fileparts(rdatafile);
tp@0:        for ii = 1:nRbatches
tp@0:            rdatafile = [rdatafilestem,'_',int2str(ii),'.mat'];
tp@0:            eval(['load ',rdatafile])
tp@0:            if ii == 1
tp@0:                allvisplanesfromr = visplanesfromr;
tp@0:            else
tp@0:               allvisplanesfromr = [allvisplanesfromr visplanesfromr]; 
tp@0:            end
tp@0:        end
tp@0:        visplanesfromr = allvisplanesfromr;
tp@0:        clear allvisplanesfromr;
tp@0:        rdatafile = rdatafilestem;
tp@0:         
tp@0:     end
tp@0:     if size(visplanesfromr,2) == 1
tp@0:         if ~any(visplanesfromr)
tp@0:             disp(' ')
tp@0:             disp('WARNING!!! The receiver can not see a single plane!')
tp@0:             disp('           Check if this is correct! ')
tp@0:             disp('           A common cause is that the receiver is placed very close to a plane,')
tp@0:             disp('           but behind it. (CR => continue calculations)')
tp@0:             disp(' ')
tp@0:             pause
tp@0:         end                
tp@0:     else
tp@0:         iv = find(sum(visplanesfromr) == 0);
tp@0:         if ~isempty(iv)
tp@0:             disp(' ')
tp@0:             disp('WARNING!!! Some of the receiver(s) can not see a single plane!')
tp@0:             disp('           These are the receiver numbers:')
tp@0:            disp(['           ',int2str(iv(:).')])
tp@0:             disp('           Check if this is correct! ')
tp@0:             disp('           A common cause is that a receiver is placed very close to a plane,')
tp@0:             disp('           but behind it. (CR => continue calculations)')
tp@0:             disp(' ')
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------------------------------------
tp@0: % The edgeseesedge test is run separately, for the cases where difforder >= 2. 
tp@0: % The data is stored in the eddata file again.
tp@0: 
tp@0: if difforder >= 2 & calcpaths == 1
tp@0:     if exist('eddata2file') ~= 1
tp@0: 	    desiredfile = [Filepath,Filestem,'_eddata'];
tp@0:     	if SHOWTEXT >= 2
tp@0: 		    disp(['   Adding edge-to-edge visibility to the eddatafile: ',eddatafile])
tp@0:     	end
tp@0:         if exist('ndiff2batches') ~= 1
tp@0:             ndiff2batches = 1;    
tp@0:         end
tp@0:     	eddatafile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches);    
tp@0: 
tp@0: 	else
tp@0: 		if SHOWTEXT >= 2
tp@0: 			disp(['   Using existing eddata2file: ',eddata2file])
tp@0: 		end
tp@0: 	end
tp@0:     
tp@0: end
tp@0: 
tp@0: 
tp@0: %---------------------------------------------------------------------------------------------
tp@0: % Make a big loop running through the source and receiver coordinates.
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: eval(['load ',eddatafile])
tp@0: clear cornerinfrontofplane planeseesplane edgeseesplane
tp@0: 
tp@0: [nplanes,slask]      = size(planecorners);
tp@0: [nedges,slask]      = size(edgecorners);
tp@0: 
tp@0: if exist('calcpaths')~=1
tp@0: 	calcpaths = 1;
tp@0: end
tp@0: if exist('calcirs')~=1
tp@0: 	calcirs = 1;
tp@0: end
tp@0: 
tp@0: % Display receiver numbers with an interval that depends on the 
tp@0: % value of SHOWTEXT.
tp@0: 
tp@0: if SHOWTEXT == 1
tp@0:     idisp = ceil(nreceivers/20);
tp@0: elseif SHOWTEXT > 1
tp@0:     idisp = 1;
tp@0: end
tp@0: 
tp@0: % We construct the string version of counters and let them be global
tp@0: % variables so we don't have to call int2str all the time
tp@0: 
tp@0: ncountersneeded = max(max([specorder difforder nplanes nedges]));
tp@0: if ncountersneeded < 10
tp@0:     JJ = setstr(32*ones(ncountersneeded,1));
tp@0: else
tp@0:     JJ = setstr(32*ones(ncountersneeded,2));
tp@0: end
tp@0: for jj=1:ncountersneeded
tp@0:    jjstr = int2str(jj);
tp@0:    JJ(jj,1:length(jjstr)) = jjstr;
tp@0: end
tp@0: [n1,n2] = size(JJ);
tp@0: JJnumbofchars = ones(n1,1);
tp@0: if n1>9
tp@0:     JJnumbofchars(10:n1) = JJnumbofchars(10:n1)+1;
tp@0:     if n1 > 99
tp@0:         JJnumbofchars(100:n1) = JJnumbofchars(100:n1)+1;    
tp@0:         if n1 > 999
tp@0:             JJnumbofchars(1000:n1) = JJnumbofchars(1000:n1)+1;    
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: NSOU = int2str(nsources);
tp@0: NREC = int2str(nreceivers);
tp@0: 
tp@0: tsetup = etime(clock,t00);
tp@0: 
tp@0: % #########################################################################
tp@0: %   
tp@0: %      Find the paths
tp@0: %   
tp@0: % #########################################################################
tp@0: 
tp@0: t00 = clock;
tp@0: tISEStree = zeros(nsources,1);
tp@0: 
tp@0: if calcpaths == 1
tp@0: 
tp@0:     if SHOWTEXT >= 1
tp@0:         disp(' ')
tp@0:         disp('####################################################################')
tp@0:         disp('#')
tp@0:         disp('#  Finding the paths')
tp@0:         disp(' ')
tp@0:     end
tp@0: 
tp@0:     % Some variables that are needed for empty combinations
tp@0: 
tp@0: 	pathtypevec = [];
tp@0: 	reflpaths = [];
tp@0: 	pathlengthvec = [];
tp@0: 	specextradata = [];
tp@0: 	edgeextradata = [];
tp@0: 	mainlistguide = [];
tp@0:     Sinsideplanenumber = [];
tp@0:     Rinsideplanenumber = [];
tp@0:     mainlistguide = [];
tp@0:     mainlistguidepattern = [];
tp@0:     directsoundrow = [];
tp@0:     allspecrows = [];
tp@0:     firstdiffrow = [];
tp@0: 	Varlist = [' pathtypevec reflpaths specextradata edgeextradata S R mainlistguide'];		
tp@0:     Varlist = [Varlist,' Sinsideplanenumber Rinsideplanenumber mainlistguide mainlistguidepattern directsoundrow allspecrows firstdiffrow'];
tp@0: 
tp@0: 	if doSbatches == 0
tp@0:         eval(['load ',sdatafile])   
tp@0: 	else
tp@0:         latestsdatafile = reftoSbatchlist(1);
tp@0:         [sdatafilepath,filestem,temp1] = fileparts(sdatafile);
tp@0:         eval(['load ',[sdatafilepath,filesep,filestem],'_',int2str(latestsdatafile),'.mat'])     
tp@0: 	end
tp@0: % % % 	% We create a souinsideplane matrix from the visplanesfroms matrix    
tp@0: % % % 	% for use in EDB1makeirs
tp@0: 	
tp@0: 	if doRbatches == 0
tp@0:         eval(['load ',rdatafile])
tp@0: 	else
tp@0:         latestrdatafile = reftoRbatchlist(1)
tp@0:         [rdatafilepath,filestem,temp1] = fileparts(rdatafile)
tp@0:         if ~isempty(rdatafilepath)
tp@0:             txtstr = ['load ',[rdatafilepath,filesep,filestem],'_',int2str(latestrdatafile),'.mat']
tp@0:         else
tp@0:             txtstr = ['load ',[filestem],'_',int2str(latestrdatafile),'.mat']            
tp@0:         end
tp@0:         eval(txtstr)     
tp@0: 	end
tp@0: % % % 	% We create a recinsideplane matrix from the visplanesfromr matrix    
tp@0: % % % 	% for use in EDB1makeirs
tp@0:     
tp@0:     for isou = soulist_for_findpaths
tp@0:     
tp@0:         t00_sou = clock;
tp@0:         ISOU = int2str(isou);
tp@0:         if SHOWTEXT >= 1
tp@0: 			disp(['Calculating for source ',ISOU,' (of ',NSOU,') '])
tp@0: 		end
tp@0:         if doSbatches == 1
tp@0:             if latestsdatafile ~= reftoSbatchlist(isou)
tp@0:                 latestsdatafile = reftoSbatchlist(isou);
tp@0:                 [sdatafilepath,filestem,temp1] = fileparts(sdatafile);
tp@0:                 eval(['load ',[sdatafilepath,filesep,filestem],'_',int2str(latestsdatafile),'.mat'])     
tp@0:             end
tp@0:             Scolnumber = isou-Sbatchlist(reftoSbatchlist(isou),1)+1
tp@0:         else
tp@0:             Scolnumber = isou;
tp@0:         end
tp@0:         S = sources(Scolnumber,:);
tp@0: 	
tp@0:         %   ###########################################################
tp@0:         %   #
tp@0:         %   #   Calculate the ISES tree for each source
tp@0:         %   #   IS = image sources, for specular reflections
tp@0:         %   #   ES = edge sources, for edge diffraction
tp@0:         %   #
tp@0:         %   ###########################################################
tp@0: 
tp@0: 	    if exist('ISEStreefile') ~= 1
tp@0:             if SHOWTEXT >= 2
tp@0:                 disp(['   Calculating ISES tree'])    
tp@0:             end
tp@0:             if difforder >= 1
tp@0:                 usedISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms(:,isou),vispartedgesfroms(:,isou),nedgesubs);
tp@0:             else
tp@0:                 usedISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms(:,isou),[],nedgesubs);
tp@0:             end
tp@0:             
tp@0:         else
tp@0:             [ISEStreefilepath,filestem,temp1] = fileparts(ISEStreefile);
tp@0:             usedISEStreefile = [[ISEStreefilepath,filesep,filestem],'_',ISOU,'_ISEStree.mat;'];   
tp@0:             if SHOWTEXT >= 2
tp@0:                 disp(['   Using existing ISEStreefile: ',ISEStreefile,'_',ISOU,'_ISEStree.mat'])    
tp@0:             end
tp@0:         end
tp@0:         if SUPPRESSFILES == 0
tp@0:             eval(['load ',usedISEStreefile])
tp@0:         else
tp@0:             POTENTIALISES = usedISEStreefile.POTENTIALISES;
tp@0:             ORIGINSFROM = usedISEStreefile.ORIGINSFROM;
tp@0:             ISCOORDS = usedISEStreefile.ISCOORDS;
tp@0:             ISESVISIBILITY = usedISEStreefile.ISESVISIBILITY;
tp@0:             IVNSPECMATRIX = usedISEStreefile. IVNSPECMATRIX;
tp@0:             lengthNspecmatrix = usedISEStreefile.lengthNspecmatrix;
tp@0:             IVNDIFFMATRIX = usedISEStreefile.IVNDIFFMATRIX;
tp@0:             lengthNdiffmatrix = usedISEStreefile.lengthNdiffmatrix;
tp@0:             singlediffcol = usedISEStreefile.singlediffcol;
tp@0:             REFLORDER = usedISEStreefile.REFLORDER;
tp@0:             startindicessinglediff = usedISEStreefile.startindicessinglediff;
tp@0:             endindicessinglediff = usedISEStreefile.endindicessinglediff;
tp@0:             ndecimaldivider = usedISEStreefile.ndecimaldivider;
tp@0:             PointertoIRcombs = usedISEStreefile.PointertoIRcombs;
tp@0:             IRoriginsfrom = usedISEStreefile.IRoriginsfrom;
tp@0:         end
tp@0:         tISEStree(isou) = etime(clock,t00_sou);
tp@0: 
tp@0:         if dobackscatter ~= 1, 
tp@0:             reclist_in_forloop = reclist_for_findpaths;
tp@0:         else
tp@0:             reclist_in_forloop = isou;            
tp@0:         end
tp@0: 
tp@0:         for irec = reclist_in_forloop
tp@0: 
tp@0:             %   ###########################################################
tp@0:             %   #
tp@0:             %   #   Find the valid paths for each source-receiver combination
tp@0:             %   #   _edpaths files are created
tp@0:             %   #
tp@0:             %   ###########################################################
tp@0: 
tp@0:             IREC = int2str(irec);
tp@0:         	if SHOWTEXT >= 1
tp@0:                 if round(irec/idisp)*idisp==irec,  
tp@0: 			        disp(['...receiver no. ',IREC,' (of ',NREC,')'])
tp@0:                 end    
tp@0:             end
tp@0:             if doRbatches == 1
tp@0:                 if latestrdatafile ~= reftoRbatchlist(irec)
tp@0:                     latestrdatafile = reftoRbatchlist(irec);
tp@0:                     [rdatafilepath,filestem,temp1] = fileparts(rdatafile);
tp@0:                     if ~isempty(rdatafilepath)
tp@0:                         txtstr = ['load ',[rdatafilepath,filesep,filestem],'_',int2str(latestrdatafile),'.mat'];
tp@0:                     else
tp@0:                         txtstr = ['load ',[filestem],'_',int2str(latestrdatafile),'.mat'];                        
tp@0:                     end
tp@0:                     eval(txtstr)     
tp@0:                 end
tp@0:                 Rcolnumber = irec-Rbatchlist(reftoRbatchlist(irec),1)+1;
tp@0:             else
tp@0:                 Rcolnumber = irec;
tp@0:             end
tp@0:             R = receivers(Rcolnumber,:);
tp@0: 
tp@0:             if exist('edpathsfile') ~= 1
tp@0: 			    desirededpathsfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0:                 if soutsidemodel(Scolnumber) == 0 & routsidemodel(Rcolnumber) == 0,        
tp@0:                     if difforder > 0
tp@0:                         usededpathsfile = EDB1findpathsISESx(eddatafile,lengthNspecmatrix,lengthNdiffmatrix,singlediffcol,startindicessinglediff,...
tp@0:                             endindicessinglediff,ndecimaldivider,PointertoIRcombs,IRoriginsfrom,...
tp@0:                             S,R,Scolnumber,Rcolnumber,directsound,specorder,difforder,...
tp@0:                             nedgesubs,visplanesfroms(:,Scolnumber),visplanesfromr(:,Rcolnumber),...
tp@0:                             vispartedgesfroms(:,Scolnumber),vispartedgesfromr(:,Rcolnumber),desirededpathsfile);
tp@0:                     else
tp@0:                         usededpathsfile = EDB1findpathsISESx(eddatafile,lengthNspecmatrix,[],[],[],...
tp@0:                             [],[],[],[],S,R,Scolnumber,Rcolnumber,directsound,specorder,difforder,...
tp@0:                             nedgesubs,visplanesfroms(:,Scolnumber),visplanesfromr(:,Rcolnumber),[],[],desirededpathsfile);
tp@0:                     end
tp@0:                 else
tp@0:                     usededpathsfile = desirededpathsfile;
tp@0:                     eval(['save ',usededpathsfile,Varlist])
tp@0:                 end
tp@0: 			else
tp@0: 				if SHOWTEXT >= 2
tp@0: 					disp(['   Using existing edpathsfile: ',edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'])
tp@0:                     desirededpathsfile = [edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0: 				end
tp@0: 			end
tp@0:             
tp@0: 	
tp@0: 		end		% of the for isou = 
tp@0: 	
tp@0: 	end		% of the for irec = 
tp@0: end
tp@0: 
tp@0: tfindpaths = etime(clock,t00);
tp@0: 
tp@0: % #########################################################################
tp@0: %   
tp@0: %      Edit paths
tp@0: %   
tp@0: % #########################################################################
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: if ~isempty(symmetricedgepairs) & exist('edpathsfile') == 1
tp@0: 
tp@0:     disp(['WARNING: You specified symmetric edge pairs, but path pruning is not'])
tp@0:     disp(['         allowed when you have specified an edpaths file'])
tp@0:     
tp@0: end
tp@0: if ~isempty(symmetricedgepairs) & exist('edpathsfile') ~= 1
tp@0: 
tp@0:     if SHOWTEXT >= 1
tp@0:         disp(' ')
tp@0:         disp('####################################################################')
tp@0:         disp('#')
tp@0:         disp('#  Pruning paths: removing redundant (symmetric) diffraction combinations')
tp@0:         disp(' ')
tp@0:     end
tp@0:         
tp@0: 	for isou = soulist_for_calcirs
tp@0:         ISOU = int2str(isou);
tp@0:         if SHOWTEXT >= 1
tp@0: 			disp(['Pruning paths for source ',ISOU,' (of ',NSOU,') '])
tp@0: 		end
tp@0: 
tp@0:         if dobackscatter ~= 1, 
tp@0:             reclist_in_forloop = reclist_for_calcirs;
tp@0:         else
tp@0:             reclist_in_forloop = isou;            
tp@0:         end
tp@0: 
tp@0:         for irec = reclist_in_forloop
tp@0:             IREC = int2str(irec);
tp@0:         	if SHOWTEXT >= 1
tp@0:                 if round(irec/idisp)*idisp==irec,  
tp@0: 			        disp(['...receiver ',IREC,' (of ',NREC,')'])
tp@0:                 end    
tp@0:             end
tp@0: 
tp@0: 			edpathsfiletoedit = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0: 
tp@0:             EDB1editpathlist(edpathsfiletoedit,[],symmetricedgepairs,[]);    
tp@0:     
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: tprunepaths = etime(clock,t00);
tp@0: 
tp@0: % #########################################################################
tp@0: %   
tp@0: %      Construct IRs
tp@0: %   
tp@0: % #########################################################################
tp@0: 
tp@0: t00 = clock;
tp@0: 
tp@0: approxplanemidpoints = 0.5*(maxvals+minvals);
tp@0: 
tp@0: if calcirs == 1
tp@0: 
tp@0:     if SHOWTEXT >= 1
tp@0:         disp(' ')
tp@0:         disp('####################################################################')
tp@0:         disp('#')
tp@0:         disp('#  Constructing IRs')
tp@0:         disp(' ')
tp@0:     end
tp@0:     
tp@0:     % If the parameter guiderowstouse wasn't specified in the setup file
tp@0:     % we give it the default value.
tp@0:     
tp@0:     if exist('guiderowstouse') ~= 1
tp@0:         guiderowstouse = [];    
tp@0:     end
tp@0: 
tp@0:     % Some variables that are needed for empty combinations
tp@0:     irgeom = [];
tp@0:     irtot = [];
tp@0:     irdirect = [];
tp@0:     irdiff = [];
tp@0:     Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR elemsize nedgesubs'];
tp@0:     Varlistaccum = [' irdirectaccum irgeomaccum irdiffaccum irtotaccum FSAMP Rstart CAIR elemsize nedgesubs'];
tp@0: 
tp@0:     reclist_in_forloop = reclist_for_calcirs;
tp@0: 
tp@0: 	for irec = reclist_in_forloop
tp@0: 		IREC = int2str(irec);
tp@0: 		if SHOWTEXT >= 1
tp@0: 			disp(['...receiver ',IREC,' (of ',NREC,')'])
tp@0: 		end
tp@0: 		
tp@0: 		irdirectaccum = [];
tp@0:         irgeomaccum = [];
tp@0:         irtotaccum = [];
tp@0:         irdiffaccum = [];
tp@0: 
tp@0: 		if dobackscatter == 1, 
tp@0: 			soulist_for_calcirs = irec;
tp@0: 		end
tp@0: 		
tp@0: 		for isou = soulist_for_calcirs
tp@0:     	    ISOU = int2str(isou);
tp@0:         	if SHOWTEXT >= 1
tp@0:                 if round(isou/idisp)*idisp==isou,  
tp@0: 					disp(['Constructing IRs for source ',ISOU,' (of ',NSOU,') '])
tp@0: 				end
tp@0: 			end
tp@0: 			
tp@0: 			if doaddsources == 1
tp@0:            		desiredirfile = [Filepath,Filestem,'_',IREC,'_ir.mat'];
tp@0:            	else
tp@0:            		desiredirfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_ir.mat'];
tp@0: 			end
tp@0: 			
tp@0:             if exist('edpathsfile') ~= 1
tp@0: 			    usededpathsfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0: 			else
tp@0: 				if SHOWTEXT >= 2
tp@0:                     if round(irec/idisp)*idisp==irec,  
tp@0: 					    disp(['   Using existing edpathsfile: ',edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'])
tp@0:                     end
tp@0: 				end
tp@0:                 usededpathsfile = [edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0: 			end
tp@0:            
tp@0:            if difforder >= 2
tp@0:                 
tp@0:       		    edirfile = EDB1makeirs(usededpathsfile,specorder,...
tp@0:                     Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0:                     edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,...
tp@0:                     reftoshortlistE,re1sho,re2sho,thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,directsound,saveindividualdiffirs);
tp@0:             else
tp@0: 
tp@0:       		    edirfile = EDB1makeirs(usededpathsfile,specorder,...
tp@0:                     Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0:                     edgelengthvec,planeeqs,[],reflfactors,closwedangvec,planesatedge,elemsize,[],[],[],[],[],[],[],[],[],desiredirfile,guiderowstouse,directsound,saveindividualdiffirs);
tp@0:                 
tp@0:             end
tp@0:  
tp@0: 	        if doaddsources == 1
tp@0: 				eval(['load ',edirfile])
tp@0: 			
tp@0: 				nirnew = length(irtot);
tp@0: 				nirold = length(irtotaccum);
tp@0: 				if nirnew > nirold
tp@0: 					irtotaccum = [irtotaccum;zeros(nirnew-nirold,1)]; 
tp@0: 					irdirectaccum = [irdirectaccum;zeros(nirnew-nirold,1)]; 
tp@0: 					irgeomaccum = [irgeomaccum;zeros(nirnew-nirold,1)]; 
tp@0: 					irdiffaccum = [irdiffaccum;zeros(nirnew-nirold,size(irdiff,2))];             
tp@0: 				end
tp@0: 				irtotaccum(1:nirnew) = irtotaccum(1:nirnew) + irtot;
tp@0: 				irdirectaccum(1:nirnew) = irdirectaccum(1:nirnew) + irdirect;
tp@0: 				irgeomaccum(1:nirnew) = irgeomaccum(1:nirnew) + irgeom;
tp@0: 				irdiffaccum(1:nirnew,:) = irdiffaccum(1:nirnew,:) + irdiff;
tp@0: 			end
tp@0: 			
tp@0:         end		% of the for isou = 1:nsou
tp@0:         
tp@0:         if doaddsources == 1
tp@0: 	        eval(['save ',edirfile,Varlistaccum])
tp@0:     	end    
tp@0: 	
tp@0: 	end		% of the for irec = 1:nrec
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------
tp@0: % Save the timing data 
tp@0: 
tp@0: tmakeirs = etime(clock,t00);
tp@0: 
tp@0: logfilename = [Filepath,Filestem,'_log.mat'];
tp@0: Varlist = [' tcadgeo tedgeo tsgeo trgeo tsetup tISEStree tfindpaths tmakeirs tprunepaths nsources nreceivers'];
tp@0: eval(['save ',logfilename,Varlist])