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