Mercurial > hg > human-echolocation

annotate private/EDB1main.m @ 18:2d5f50205527 jabuilder_int tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Escape the trailing backslash as well
author Chris Cannam
date Tue, 30 Sep 2014 16:23:00 +0100
parents d9262cdbfb38
children
rev   line source
tp@0 1 function EDB1main(EDsetupfile)
tp@0 2 % EDB1main - Calculates the specular and edge diffraction IRs and saves them in a file.
tp@0 3 % Calculates the specular and edge diffraction IRs and saves them in a file or
tp@0 4 % a series of files, if there are many sources and receivers.The
tp@0 5 % calculation goes through three stages:
tp@0 6 % 1. Geometrical precalculations: define edges and determine visibility
tp@0 7 % between edges, planes, sources and receivers.
tp@0 8 % 2. Find the valid sound paths and save a list of them
tp@0 9 % 3. Go through the list and construct IRs
tp@0 10 %
tp@0 11 % This version, EDB1main, includes the option of dividing the
tp@0 12 % calculations into batches. When there are many receivers, the geometrical
tp@0 13 % precalculations might overload the available memory since all receivers
tp@0 14 % are treated simulataneously for each source.
tp@0 15 %
tp@0 16 % Input parameters:
tp@0 17 % EDsetupfile (optional) This file contains all parameter settings for the calculations.
tp@0 18 % If no filename is given, a file open window will appear and the
tp@0 19 % desired EDsetupfile can be chosen.
tp@0 20 % The file could either be a .m file (the "EDsetupfile")
tp@0 21 % or a .mat file (the "setupmatfile").
tp@0 22 % Input parameters in the EDsetupfile:
tp@0 23 % FSAMP, CAIR, RHOAIR, SHOWTEXT
tp@0 24 % (global) The sampling frequency, speed of sound, density of air
tp@0 25 % and the control of how much text should
tp@0 26 % be plotted (SHOWTEXT could be 0,1,2,3,4).
tp@0 27 % SUPPRESSFILES (global, optional) If this variable is introduced in the
tp@0 28 % setupfile as a global variable, with the value 1, then the files
tp@0 29 % 'xxx_ISEStree.mat' will not be generated. Instead values are passed
tp@0 30 % directly.
tp@0 31 % Filepath The path to the folder where all output data files will be stored.
tp@0 32 % Filestem The first part of the output data files will be given this name
tp@0 33 % CADfile The name of the .CAD file where the geometry is given, including path
tp@0 34 % open_or_closed_model Specifies whether the geometrical model is open or closed:
tp@0 35 % 'o' Open model
tp@0 36 % 'c' Closed model
tp@0 37 % int_or_ext_model Specifies whether the geometrical model is interior or exterior:
tp@0 38 % 'i' Interior model (NB! An interior model must be closed,
tp@0 39 % otherwise it should be called an exterior model)
tp@0 40 % 'e' Exterior model
tp@0 41 % EDcalcmethod Specifies the edge diffraction method. Must have one of the following values:
tp@0 42 % 'v' Vanderkooy
tp@0 43 % 'n' Svensson et al
tp@0 44 % 'k' Kirchhoff approximation NB! The Kirchhoff
tp@0 45 % approximation gives only first-order diffraction.
tp@0 46 % directsound 0 or 1
tp@0 47 % specorder The maximum order of specular reflections. (To get all possible combinations
tp@0 48 % specorder and difforder should be given the same value).
tp@0 49 % difforder The maximum order of diffraction components. (difforder can never be higher
tp@0 50 % than specorder, i.e., if you want difforder = 2, you must set specorder to 2 or higher)
tp@0 51 % elemsize A vector of size [1,difforder] which specifies the size of the edge subdivision
tp@0 52 % for edge diffraction calculations. A value of 0.2 means that the edge elements
tp@0 53 % are 5 times (1/0.2 =5) the wavelength at FSAMP.
tp@0 54 % A higher value gives more accurate edge diffraction but is also more time consuming.
tp@0 55 % Larger and larger elements could be employed for higher-order diffraction.
tp@0 56 % Suggested values for high accuracy: elemsize = [1 0.5 0.25 0.125 ....].
tp@0 57 % nedgesubs All edges are subdivided very coarsely into edge segments, for visibility checks.
tp@0 58 % The value of nedgesubs states how many such segments each edge will be divided into.
tp@0 59 % The minimum value is 2. NB! All edges will be divided into equally many segments
tp@0 60 % in order to employ efficient Matlab vector-based processing.
tp@0 61 % sources A vector of the source coordinates, size [nsources,3], or source names/numbers
tp@0 62 % size [nsources,1] or [nsources,2] (if they should be taken from the CADfile.)
tp@0 63 % receivers A vector of the receiver coordinates or receiver numbers (if they should be taken
tp@0 64 % from the CADfile. Same syntax as for sources.
tp@0 65 % nSbatches,nRbatches
tp@0 66 % (optional) An integer >= 1 which defines how many chunks the number of sources/receivers
tp@0 67 % will be divided into. Default = 1. nSbatches/nRbatches = 1 means that all sources/receivers are
tp@0 68 % treated simultaneously. nSbatches/nRbatches = 2 means that the list of sources/receivers
tp@0 69 % is divided into two chunks etc. It is fastest with nSbatches/nRbatches = 1 but it also
tp@0 70 % requires the largest amount of memory. (For large problems, there might even be no time
tp@0 71 % penalty in choosing nRbatches = the number of receivers!)
tp@0 72 % soulist_for_findpaths/reclist_for_findpaths
tp@0 73 % (optional) A list of integers that specifies which of the sources/receivers that
tp@0 74 % should be run for the "find paths" stage. Default = [1:nsources]/[1:nreceivers].
tp@0 75 % This is useful for huge problems where one might want to stop (with CTRL-C) and restart calculations.
tp@0 76 % soulist_for_calcirs/reclist_for_calcirs
tp@0 77 % (optional) A list of integers that specifies which of the sources/receivers that
tp@0 78 % should be run for the "calc irs" stage. Default = [1:nsources]/[1:nreceivers].
tp@0 79 % This is useful for huge problems where one might want to stop (with CTRL-C) and restart calculations.
tp@0 80 % firstskipcorner (optional) All edges including at least one corner with this number or
tp@0 81 % higher will be excluded from the calculations. Default: 1e6
tp@0 82 % Rstart (optional) The distance that should correspond to the time zero for
tp@0 83 % all IRs. Default = 0. NB! With Rstart = 0, all IRs will include the initial time
tp@0 84 % delay that is caused by the propagation time from the source to the receiver. For
tp@0 85 % long distances and high FSAMP, this can lead to many initial samples that are zero.
tp@0 86 % If you know that all source-receiver combinations have a certain minimum
tp@0 87 % distance, this minimum distance (or a value slightly below) can be given as Rstart.
tp@0 88 % If Rstart is given a value which is smaller than the shortest distance encountered,
tp@0 89 % an error will occur.
tp@0 90 %
tp@0 91 % Optional input parameters in the EDsetupfile:
tp@0 92 % cadgeofile If this parameter is specified, this file will be read rather than the CADfile.
tp@0 93 % Saves some time for big geometries.
tp@0 94 % eddatafile If this parameter is specified, the function EDB1edgeo is not called. Saves
tp@0 95 % some time.
tp@0 96 % srdatafile If this parameter is specified, the function EDB1srgeo is not called. Saves
tp@0 97 % some time.
tp@0 98 % restartnumber If this parameter is specified, the calculations will start at receiver number
tp@0 99 % restartnumber, rather than number 1. This is useful if a calculation for many
tp@0 100 % receivers can not be completed in one run.
tp@0 101 % stopnumber If this parameter is specified, the calculations will stop for receiver number
tp@0 102 % stopnumber.
tp@0 103 % saveindividualdiffirs A vector with two values, 0 or 1.
tp@0 104 % If the first value is given the value 1, the diffraction
tp@0 105 % IR will have one column for each order. Default is that
tp@0 106 % all orders are summed into a single irdiff.
tp@0 107 % If the second value is given the value 1, all individual diffraction irs
tp@0 108 % will be saved on individual lines in a large matrix called 'alldiffirs'.
tp@0 109 % This matrix has only single precision.
tp@0 110 % Default value: [0 0]
tp@0 111 % doaddsources 0 or 1, (default value 0) indicating whether:
tp@0 112 % 0: every individual S-to-R IR should be saved
tp@0 113 % 1: a single IR is saved for each R, which is
tp@0 114 % constructed by first computing all individual
tp@0 115 % S-to-R IRs, saving them in a temp file, and then
tp@0 116 % add IRs from all sources together. This way a
tp@0 117 % loudspeaker element can be simulated by a
tp@0 118 % distribution of point sources.
tp@0 119 % dobackscatter 0 or 1, (default value 0), indicating whether:
tp@0 120 % 0: every individual S-to-R IR should be saved
tp@0 121 % 1: only S1 to R1, S2 to R2, S3 to R3 combos
tp@0 122 % are calculated.
tp@0 123 %
tp@0 124 % Output parameters:
tp@0 125 % The following IRs are calculated and saved in the output file(s). Any of them can be
tp@0 126 % excluded from the calculations by setting the parameters in the setup file.
tp@0 127 %
tp@0 128 % irtot the total impulse response
tp@0 129 % irdirect the direct sound
tp@0 130 % irgeom the purely geometrical acoustics components, i.e. specular reflections.
tp@0 131 % irdiff the diffracted components. This will innclude all combinations of specular and
tp@0 132 % diffraction, and first-order as well as higher-order combinations.
tp@0 133 %
tp@0 134 % Note that these variables, together with those in the EDsetupfile are saved in a results file
tp@0 135 % which gets the name in Filestem (see below) + '_SS_RR_ir.mat' (if addsources = 0) or
tp@0 136 % '_RR_ir.mat' (if addsources = 1).
tp@0 137 %
tp@0 138 % Note also that the impulse responses are scaled such that the direct sound gets an amplitude of
tp@0 139 % 1/distance. This means that the impulse response represents a system with the sound pressure as
tp@0 140 % output signal and an input signal = RHOAIR*Volumeacc(t)/(4*pi) where Volumeacc(t) is
tp@0 141 % the volume acceleration of the monopole source.
tp@0 142 %
tp@0 143 % A logfile containing the computation time for the three stages is
tp@0 144 % created, with the name Filestem + '_log.mat'
tp@0 145 %
tp@0 146 % Uses functions EDB1readcad EDB1readsetup EDB1edgeo
tp@0 147 % EDB1Sgeo EDB1Rgeo EDB1findISEStree EDB1findpathsISES EDB1makeirs
tp@0 148 %
tp@0 149 % ----------------------------------------------------------------------------------------------
tp@0 150 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 151 %
tp@0 152 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 153 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 154 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 155 %
tp@0 156 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 157 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 158 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 159 %
tp@0 160 % You should have received a copy of the GNU General Public License along with the
tp@0 161 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 162 % ----------------------------------------------------------------------------------------------
tp@0 163 % Peter Svensson (svensson@iet.ntnu.no) 20100211
tp@0 164 %
tp@0 165 % EDB1main(EDsetupfile);
tp@0 166
tp@0 167 global FSAMP CAIR RHOAIR SHOWTEXT JJ JJnumbofchars SUPPRESSFILES
tp@0 168 global POTENTIALISES ISCOORDS IVNDIFFMATRIX
tp@0 169 global IVNSPECMATRIX ORIGINSFROM ISESVISIBILITY REFLORDER
tp@0 170
tp@0 171 global IRDIFFVEC
tp@0 172
tp@0 173 FSAMP = []; CAIR = []; RHOAIR = []; SHOWTEXT = [];
tp@0 174
tp@0 175 %--------------------------------------------------------------
tp@0 176 % First run the setup m file, just to get access to the CAD-file name
tp@0 177 % so that the CAD-file can be converted into a cadgeofile.
tp@0 178 % (This conversion must be done before the setup m-file is converted
tp@0 179 % to a setup mat-file since the source and receiver coordinates should
tp@0 180 % be specified in the setupmatfile, and the sources and receivers might
tp@0 181 % have to be read from the CAD-file).
tp@0 182 %
tp@0 183 % First we make a little check of a common mistake: check if the paths
tp@0 184 % seem to be for the correct computer type.
tp@0 185 % The EDsetupfile-path does not need to be checked if a file window is
tp@0 186 % presented. The CAD-file path always needs to be checked.
tp@0 187
tp@0 188 % We check if the type of setupfile:
tp@0 189 % .m setupfiletype = 1;
tp@0 190 % .mat setupfiletype = 2;
tp@0 191
tp@0 192 setupfiletype = 1;
tp@0 193
tp@0 194 if nargin == 0
tp@0 195 [EDsetupfile,EDsetupfilepath] = uigetfile('*.m','Please select the EDsetupfile');
tp@0 196 [temp1,EDsetupfile,temp2] = fileparts([EDsetupfilepath,EDsetupfile]);
tp@0 197
tp@0 198 EDsetupfile = [EDsetupfilepath,EDsetupfile];
tp@0 199
tp@0 200 if ~isstr(EDsetupfile) | isempty(EDsetupfile)
tp@0 201 return
tp@0 202 end
tp@0 203 else
tp@0 204 [EDsetupfilepath,tempfilestem,fileext] = fileparts(EDsetupfile);
tp@0 205 EDsetupfilepath = [EDsetupfilepath,filesep];
tp@0 206 if length(fileext) == 4
tp@0 207 if fileext == '.mat'
tp@0 208 setupfiletype = 2;
tp@0 209 end
tp@0 210 end
tp@0 211 end
tp@0 212
luis@9 213 % need to read setup file manually --luisf.
tp@0 214
luis@9 215 % if setupfiletype == 1
luis@9 216 % comptype = computer;
luis@9 217 % if length(comptype) == 3
luis@9 218 % comptype = [comptype,'P'];
luis@9 219 % end
luis@9 220 % if comptype(1:4) == 'MACI',
luis@9 221 % [temp1,tempfilestem,temp2] = fileparts(EDsetupfile);
luis@9 222 % eval([tempfilestem])
luis@9 223 % else
luis@9 224 % eval(['run ''',EDsetupfile,''''])
luis@9 225 % end
luis@9 226 % else
luis@9 227 % eval(['load ''',EDsetupfile,''''])
luis@9 228 % end
luis@9 229
luis@9 230 eval(char(textread(EDsetupfile, '%s', 'whitespace', '')));
luis@9 231 % ~luisf
tp@0 232
tp@0 233 if SHOWTEXT >= 1
tp@0 234 disp(' ')
tp@0 235 disp('####################################################################')
tp@0 236 disp('#')
tp@0 237 disp('# EDB1mainISESx, version 11 February 2010')
tp@0 238 disp('#')
tp@0 239 end
tp@0 240
tp@0 241 if SHOWTEXT >= 2
tp@0 242 disp('#')
tp@0 243 disp('# This version does the following:')
tp@0 244 disp('# * Caclulates impulse responses that include specular reflections up to any order and')
tp@0 245 disp('# multiple diffraction up to order 6. Combinations are also taken into account, such as')
tp@0 246 disp('# specular-diffraction-specular but higher-order diffractions must be')
tp@0 247 disp('# in a single sequence, possibly with specular reflections before and after.')
tp@0 248 disp('# There will be memory problems for higher-order specular reflections due to a')
tp@0 249 disp('# memory-hungry vectorized version of the image source method.')
tp@0 250 disp('# * Combinations with two diffractions that have specular reflections inbetween')
tp@0 251 disp('# are handled, but not with any specular reflections after the last diffraction.')
tp@0 252 disp('# * Partly visible edges are handled in first-order diffraction')
tp@0 253 disp('# calculations. The accuracy depends on how high value of')
tp@0 254 disp('# nedgesubs that is chosen. Edges that are split into two or more subedges are handled properly.')
tp@0 255 disp('# * Specular reflections in the IR are handled by splitting a pulse between two consecutive time slots.')
tp@0 256 disp('# This gives quite a severe low-pass filter roll-off, but zero phase-error.')
tp@0 257 disp('# Choose a high over-sampling ratio to minimize this magnitude error.')
tp@0 258 disp('# In addition, the low-pass filter effect can be avoided altogether for the direct sound by setting')
tp@0 259 disp('# Rstart to exactly the source-receiver distance. Thereby the direct sound pulse will occur exactly in the')
tp@0 260 disp('# first time sample. For more tweaking, the three values Rstart and CAIR and FSAMP should be possible')
tp@0 261 disp('# to adjust so that one, two or three geometrical pulses arrive exactly at sample times.')
tp@0 262 disp('# * An accurate integration technique is used for first-order diffraction so that receiver positions close')
tp@0 263 disp('# to zone boundaries are handled without problem. Exactly at the zone boundary the specular reflection')
tp@0 264 disp('# or direct sound gets half the amplitude (See below though).')
tp@0 265 disp('#')
tp@0 266 disp('# Limitations with this version:')
tp@0 267 disp('# 1. Some, but not all, possible combinations where two or more diffracted')
tp@0 268 disp('# reflections have specular reflections in-between are calculated. ')
tp@0 269 disp('# This requires a bit of work to complete.')
tp@0 270 disp('# 3. In-plane visibility of edges is only partly checked, i.e., if a plane has indents')
tp@0 271 disp('# so that an edge can not see all other edges of its own plane.')
tp@0 272 disp('# 4. Part-visibility of higher-order diffraction is not checked. These components are treated')
tp@0 273 disp('# as fully visible or totally invisible. However, if the first or last edge in the')
tp@0 274 disp('# sequence is partly visible, this will be used.')
tp@0 275 disp('# 5. Non-rigid specular reflections or edge diffraction for non-rigid planes are not implemented. ')
tp@0 276 disp('# Pressure-release surfaces would be quite easy to implement.')
tp@0 277 disp('# Other surface impedances could be implemented since the hitpoints for all specular')
tp@0 278 disp('# reflections are stored in the edpaths output file.')
tp@0 279 disp('# 6. Other receivers than omnidirectional microphones are not implemented. HRTFs or directional microphones')
tp@0 280 disp('# could be implemented since the hitpoints for all specular reflections are stored in the edpaths')
tp@0 281 disp('# so the incidence angles are easy to calculate. Reflections that end with a diffraction need some')
tp@0 282 disp('# more consideration.')
tp@0 283 disp('# 7. The special cases where the receiver is exactly at a zone boundary is handled correctly only')
tp@0 284 disp('# for first-order diffraction. This takes a bit of work to implement for higher-order diffraction.')
tp@0 285 disp('# 8. Sources and receivers can not be place directly at a surface. They must be raised by 1 mm or so.')
tp@0 286 disp('#')
tp@0 287 disp('# This version could be made more efficient by:')
tp@0 288 disp('# 1. Employing ray-tracing or beam-tracing for finding the visibility of')
tp@0 289 disp('# plane to plane etc. Now, a simple plane-to-plane visibility check')
tp@0 290 disp('# is performed: are planes in front of each other? The final visibility')
tp@0 291 disp('# check is done for each image source to each receiver. An improved method would be more efficient')
tp@0 292 disp('# for geometries with lots of indents and obscuring planes, e.g., as in a city')
tp@0 293 disp('# geometry.')
tp@0 294 disp('# 2. Developing an algorithm that calculates the diffraction IR in an iterative process')
tp@0 295 disp('# instead of restarting for each new order of diffraction.')
tp@0 296 disp('# 3. The direct sound visibility test could be done earlier (in EDB1SRgeo).')
tp@0 297 disp('# 4. Calculation of Image Receiver coordinates could probably be made more efficiently.')
tp@0 298 disp('#')
tp@0 299 end
tp@0 300
tp@0 301 if SHOWTEXT >= 1
tp@0 302 disp(' ')
tp@0 303 disp([' Using the setupfile: ',EDsetupfile])
tp@0 304 disp(' ')
tp@0 305 end
tp@0 306
tp@0 307 % #########################################################################
tp@0 308 %
tp@0 309 % Geometry pre-calculations
tp@0 310 %
tp@0 311 % #########################################################################
tp@0 312
tp@0 313 %---------------------------------------------------------------------------------------------
tp@0 314 % Convert the CAD-file into a mat file - or read the cadgeofile if it is specified
tp@0 315
tp@0 316 t00 = clock;
tp@0 317
tp@0 318 if SHOWTEXT >= 1
tp@0 319 disp(' ')
tp@0 320 disp('####################################################################')
tp@0 321 disp('#')
tp@0 322 disp('# Geometry pre-calculations')
tp@0 323 disp(' ')
tp@0 324 end
tp@0 325
tp@0 326 if exist('cadgeofile') ~= 1
tp@0 327 desiredname = [Filepath,Filestem,'_cadgeo'];
tp@0 328 if SHOWTEXT >= 2
tp@0 329 disp([' Reading CAD-file: ',CADfile,' and converting it to a cadgeofile: ',desiredname])
tp@0 330 end
tp@0 331 % Edited by APO: Added support for .ac files (only geometry)
tp@0 332 [temp1,temp2,fileext] = fileparts(CADfile);
tp@0 333 if(strcmp(fileext,'.ac'))
tp@0 334 cadgeofile = EDB1readac(CADfile,desiredname);
tp@0 335 else
tp@0 336 cadgeofile = EDB1readcad(CADfile,desiredname);
tp@0 337 end
tp@0 338 else
tp@0 339 if SHOWTEXT >= 2
tp@0 340 disp([' Using existing cadgeofile: ',cadgeofile])
tp@0 341 end
tp@0 342 end
tp@0 343
tp@0 344 %--------------------------------------------------------------
tp@0 345 % Convert the EDsetupfile into a setupmatfile. If sources and receivers should
tp@0 346 % be read in the CADfile, add these to the setupmatfile.
tp@0 347 %
tp@0 348 % This is needed only if the setupfile was of type 1, that is, a .m file.
tp@0 349
tp@0 350 if setupfiletype == 1
tp@0 351 setupmatfile = EDB1readsetup(EDsetupfile,cadgeofile);
tp@0 352
tp@0 353 eval(['load ',setupmatfile])
tp@0 354 end
tp@0 355
tp@0 356 tcadgeo = etime(clock,t00);
tp@0 357
tp@0 358 %--------------------------------------------------------------
tp@0 359 % Derive the edge parameters.
tp@0 360
tp@0 361 t00 = clock;
tp@0 362
tp@0 363 if calcpaths == 1
tp@0 364
tp@0 365 if exist('eddatafile') ~= 1
tp@0 366 desiredfile = [Filepath,Filestem,'_eddata'];
tp@0 367 if SHOWTEXT >= 2
tp@0 368 disp([' Calculating edges and creating an eddatafile: ',desiredfile])
tp@0 369 end
tp@0 370 eddatafile = EDB1edgeox(cadgeofile,desiredfile,lower(open_or_closed_model(1)),lower(int_or_ext_model(1)),...
tp@0 371 specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy);
tp@0 372 else
tp@0 373 if SHOWTEXT >= 2
tp@0 374 disp([' Using existing eddatafile: ',eddatafile])
tp@0 375 end
tp@0 376 end
tp@0 377 else
tp@0 378 if exist('eddatafile') ~= 1
tp@0 379 desiredfile = [Filepath,Filestem,'_eddata'];
tp@0 380 eddatafile = desiredfile;
tp@0 381 else
tp@0 382 if SHOWTEXT >= 2
tp@0 383 disp([' Using existing eddatafile: ',eddatafile])
tp@0 384 end
tp@0 385 end
tp@0 386 end
tp@0 387
tp@0 388 tedgeo = etime(clock,t00);
tp@0 389
tp@0 390 %--------------------------------------------------------------
tp@0 391 % If batches should be used, divide the list of receivers into batches.
tp@0 392
tp@0 393 nsources = size(sources,1);
tp@0 394 nreceivers = size(receivers,1);
tp@0 395
tp@0 396 if nSbatches == 1
tp@0 397 doSbatches = 0;
tp@0 398 elseif nSbatches > 1
tp@0 399 if nSbatches > nsources
tp@0 400 error('ERROR: nSbatches can not be larger than the number of sources')
tp@0 401 end
tp@0 402 doSbatches = 1;
tp@0 403 Sbatchlist = zeros(nSbatches,2);
tp@0 404 nperbatch = ceil(nsources/nSbatches);
tp@0 405 Sbatchlist(:,2) = [1:nSbatches].'*nperbatch;
tp@0 406 Sbatchlist(:,1) = [0:nSbatches-1].'*nperbatch+1;
tp@0 407 Sbatchlist(nSbatches,2) = nsources;
tp@0 408 reftoSbatchlist = [1:nsources].';
tp@0 409 reftoSbatchlist = reftoSbatchlist(:,ones(1,nSbatches));
tp@0 410 compvec = Sbatchlist(:,1).';
tp@0 411 reftoSbatchlist = reftoSbatchlist >= compvec(ones(nsources,1),:);
tp@0 412 reftoSbatchlist = sum(reftoSbatchlist.').';
tp@0 413 end
tp@0 414
tp@0 415 if nRbatches == 1
tp@0 416 doRbatches = 0;
tp@0 417 elseif nRbatches > 1
tp@0 418 if nRbatches > nreceivers
tp@0 419 error('ERROR: nRbatches can not be larger than the number of receivers')
tp@0 420 end
tp@0 421 doRbatches = 1;
tp@0 422 Rbatchlist = zeros(nRbatches,2);
tp@0 423 nperbatch = ceil(nreceivers/nRbatches);
tp@0 424 Rbatchlist(:,2) = [1:nRbatches].'*nperbatch;
tp@0 425 Rbatchlist(:,1) = [0:nRbatches-1].'*nperbatch+1;
tp@0 426 Rbatchlist(nRbatches,2) = nreceivers;
tp@0 427 reftoRbatchlist = [1:nreceivers].';
tp@0 428 reftoRbatchlist = reftoRbatchlist(:,ones(1,nRbatches));
tp@0 429 compvec = Rbatchlist(:,1).';
tp@0 430 reftoRbatchlist = reftoRbatchlist >= compvec(ones(nreceivers,1),:);
tp@0 431 reftoRbatchlist = sum(reftoRbatchlist.').';
tp@0 432 end
tp@0 433
tp@0 434 %---------------------------------------------------------------------------------------------
tp@0 435 % Find out the S and R related parameters, make obstruction checks etc by calling
tp@0 436 % EDB1Sgeo/EDB1Rgeo or load existing sdata/rdata files.
tp@0 437
tp@0 438 t00 = clock;
tp@0 439
tp@0 440 if calcpaths == 1
tp@0 441 if exist('sdatafile') ~= 1
tp@0 442 if doSbatches == 0
tp@0 443 desiredname = [Filepath,Filestem,'_sdata'];
tp@0 444 if SHOWTEXT >= 2
tp@0 445 disp([' Calculating S parameters and creating an sdatafile: ',desiredname])
tp@0 446 end
tp@0 447
tp@0 448 sdatafile = EDB1SorRgeo(eddatafile,desiredname,sources,'S',difforder,nedgesubs);
tp@0 449 else
tp@0 450 error(['ERROR: S-batches not implemented yet'])
tp@0 451 for ii = 1:nSbatches
tp@0 452 desiredname = [Filepath,Filestem,'_sdata'];
tp@0 453 if SHOWTEXT >= 2
tp@0 454 disp([' Calculating S parameters and creating an sdatafile: ',desiredname,'_',int2str(ii)])
tp@0 455 end
tp@0 456 sdatafile = EDB1SorRgeo(eddatafile,[desiredname,'_',int2str(ii)],sources(Sbatchlist(ii,1):Sbatchlist(ii,2),:),'S',difforder,nedgesubs);
tp@0 457 sdatafile = desiredname; % To trim off the _1 or _2 or...
tp@0 458 end
tp@0 459 end
tp@0 460 else
tp@0 461 if SHOWTEXT >= 2
tp@0 462 disp([' Using existing sdatafile: ',sdatafile])
tp@0 463 end
tp@0 464 end
tp@0 465 else
tp@0 466 if exist('sdatafile') ~= 1
tp@0 467 desiredname = [Filepath,Filestem,'_sdata'];
tp@0 468 sdatafile = desiredname;
tp@0 469 else
tp@0 470 if SHOWTEXT >= 2
tp@0 471 disp([' Using existing sdatafile: ',sdatafile])
tp@0 472 end
tp@0 473 end
tp@0 474 end
tp@0 475
tp@0 476 tsgeo = etime(clock,t00);
tp@0 477
tp@0 478 %------------------------------------
tp@0 479 % We check that there is not a small geometrical mistake which makes the
tp@0 480 % source be behind all planes - but only if the user has set SHOWTEXT.
tp@0 481
tp@0 482 if SHOWTEXT >= 2
tp@0 483 eval(['load ',sdatafile])
tp@0 484 if size(visplanesfroms,2) == 1
tp@0 485 if ~any(visplanesfroms)
tp@0 486 disp(' ')
tp@0 487 disp('WARNING!!! The source(s) can not see a single plane!')
tp@0 488 disp(' Check if this is correct! ')
tp@0 489 disp(' A common cause is that the source is placed very close to a plane,')
tp@0 490 disp(' but behind it. (CR => continue calculations)')
tp@0 491 disp(' ')
tp@0 492 pause
tp@0 493 end
tp@0 494 else
tp@0 495 iv = find(sum(visplanesfroms) == 0);
tp@0 496 if ~isempty(iv)
tp@0 497 disp(' ')
tp@0 498 disp('WARNING!!! Some of the source(s) can not see a single plane!')
tp@0 499 disp(' These are the sources number:')
tp@0 500 disp([' ',int2str(iv(:).')])
tp@0 501 disp(' Check if this is correct! ')
tp@0 502 disp(' A common cause is that the source is placed very close to a plane,')
tp@0 503 disp(' but behind it. (CR => continue calculations)')
tp@0 504 disp(' ')
tp@0 505 end
tp@0 506 end
tp@0 507 end
tp@0 508
tp@0 509 t00 = clock;
tp@0 510
tp@0 511 if calcpaths == 1
tp@0 512 if exist('rdatafile') ~= 1
tp@0 513 if doRbatches == 0
tp@0 514 desiredname = [Filepath,Filestem,'_rdata'];
tp@0 515 if SHOWTEXT >= 2
tp@0 516 disp([' Calculating R parameters and creating an rdatafile: ',desiredname])
tp@0 517 end
tp@0 518 rdatafile = EDB1SorRgeo(eddatafile,desiredname,receivers,'R',difforder,nedgesubs);
tp@0 519 else
tp@0 520 for ii = 1:nRbatches
tp@0 521 desiredname = [Filepath,Filestem,'_rdata'];
tp@0 522 if SHOWTEXT >= 2
tp@0 523 disp([' Calculating R parameters and creating an rdatafile: ',desiredname,'_',int2str(ii)])
tp@0 524 end
tp@0 525
tp@0 526 rdatafile = EDB1SorRgeo(eddatafile,[desiredname,'_',int2str(ii)],receivers(Rbatchlist(ii,1):Rbatchlist(ii,2),:),'R',difforder,nedgesubs);
tp@0 527 rdatafile = desiredname; % To trim off the _1 or _2 or...
tp@0 528 end
tp@0 529 end
tp@0 530 else
tp@0 531 if SHOWTEXT >= 2
tp@0 532 disp([' Using existing rdatafile: ',rdatafile])
tp@0 533 end
tp@0 534 end
tp@0 535 else
tp@0 536 if exist('rdatafile') ~= 1
tp@0 537 desiredname = [Filepath,Filestem,'_rdata'];
tp@0 538 rdatafile = desiredname;
tp@0 539 else
tp@0 540 if SHOWTEXT >= 2
tp@0 541 disp([' Using existing rdatafile: ',rdatafile])
tp@0 542 end
tp@0 543 end
tp@0 544 end
tp@0 545
tp@0 546 trgeo = etime(clock,t00);
tp@0 547
tp@0 548 %------------------------------------
tp@0 549 % We check that there is not a small geometrical mistake which makes the
tp@0 550 % receiver(s) be behind all planes - but only if the user has set SHOWTEXT.
tp@0 551
tp@0 552 if SHOWTEXT >= 2
tp@0 553 if exist('nRbatches') == 0
tp@0 554 nRbatches = 1;
tp@0 555 end
tp@0 556 if nRbatches == 1
tp@0 557 eval(['load ',rdatafile])
tp@0 558 else
tp@0 559 [temp1,rdatafilestem,temp2] = fileparts(rdatafile);
tp@0 560 for ii = 1:nRbatches
tp@0 561 rdatafile = [rdatafilestem,'_',int2str(ii),'.mat'];
tp@0 562 eval(['load ',rdatafile])
tp@0 563 if ii == 1
tp@0 564 allvisplanesfromr = visplanesfromr;
tp@0 565 else
tp@0 566 allvisplanesfromr = [allvisplanesfromr visplanesfromr];
tp@0 567 end
tp@0 568 end
tp@0 569 visplanesfromr = allvisplanesfromr;
tp@0 570 clear allvisplanesfromr;
tp@0 571 rdatafile = rdatafilestem;
tp@0 572
tp@0 573 end
tp@0 574 if size(visplanesfromr,2) == 1
tp@0 575 if ~any(visplanesfromr)
tp@0 576 disp(' ')
tp@0 577 disp('WARNING!!! The receiver can not see a single plane!')
tp@0 578 disp(' Check if this is correct! ')
tp@0 579 disp(' A common cause is that the receiver is placed very close to a plane,')
tp@0 580 disp(' but behind it. (CR => continue calculations)')
tp@0 581 disp(' ')
tp@0 582 pause
tp@0 583 end
tp@0 584 else
tp@0 585 iv = find(sum(visplanesfromr) == 0);
tp@0 586 if ~isempty(iv)
tp@0 587 disp(' ')
tp@0 588 disp('WARNING!!! Some of the receiver(s) can not see a single plane!')
tp@0 589 disp(' These are the receiver numbers:')
tp@0 590 disp([' ',int2str(iv(:).')])
tp@0 591 disp(' Check if this is correct! ')
tp@0 592 disp(' A common cause is that a receiver is placed very close to a plane,')
tp@0 593 disp(' but behind it. (CR => continue calculations)')
tp@0 594 disp(' ')
tp@0 595 end
tp@0 596 end
tp@0 597 end
tp@0 598
tp@0 599 %---------------------------------------------------------------------------------------------
tp@0 600 % The edgeseesedge test is run separately, for the cases where difforder >= 2.
tp@0 601 % The data is stored in the eddata file again.
tp@0 602
tp@0 603 if difforder >= 2 & calcpaths == 1
tp@0 604 if exist('eddata2file') ~= 1
tp@0 605 desiredfile = [Filepath,Filestem,'_eddata'];
tp@0 606 if SHOWTEXT >= 2
tp@0 607 disp([' Adding edge-to-edge visibility to the eddatafile: ',eddatafile])
tp@0 608 end
tp@0 609 if exist('ndiff2batches') ~= 1
tp@0 610 ndiff2batches = 1;
tp@0 611 end
tp@0 612 eddatafile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches);
tp@0 613
tp@0 614 else
tp@0 615 if SHOWTEXT >= 2
tp@0 616 disp([' Using existing eddata2file: ',eddata2file])
tp@0 617 end
tp@0 618 end
tp@0 619
tp@0 620 end
tp@0 621
tp@0 622
tp@0 623 %---------------------------------------------------------------------------------------------
tp@0 624 % Make a big loop running through the source and receiver coordinates.
tp@0 625
tp@0 626 t00 = clock;
tp@0 627
tp@0 628 eval(['load ',eddatafile])
tp@0 629 clear cornerinfrontofplane planeseesplane edgeseesplane
tp@0 630
tp@0 631 [nplanes,slask] = size(planecorners);
tp@0 632 [nedges,slask] = size(edgecorners);
tp@0 633
tp@0 634 if exist('calcpaths')~=1
tp@0 635 calcpaths = 1;
tp@0 636 end
tp@0 637 if exist('calcirs')~=1
tp@0 638 calcirs = 1;
tp@0 639 end
tp@0 640
tp@0 641 % Display receiver numbers with an interval that depends on the
tp@0 642 % value of SHOWTEXT.
tp@0 643
tp@0 644 if SHOWTEXT == 1
tp@0 645 idisp = ceil(nreceivers/20);
tp@0 646 elseif SHOWTEXT > 1
tp@0 647 idisp = 1;
tp@0 648 end
tp@0 649
tp@0 650 % We construct the string version of counters and let them be global
tp@0 651 % variables so we don't have to call int2str all the time
tp@0 652
tp@0 653 ncountersneeded = max(max([specorder difforder nplanes nedges]));
tp@0 654 if ncountersneeded < 10
tp@0 655 JJ = setstr(32*ones(ncountersneeded,1));
tp@0 656 else
tp@0 657 JJ = setstr(32*ones(ncountersneeded,2));
tp@0 658 end
tp@0 659 for jj=1:ncountersneeded
tp@0 660 jjstr = int2str(jj);
tp@0 661 JJ(jj,1:length(jjstr)) = jjstr;
tp@0 662 end
tp@0 663 [n1,n2] = size(JJ);
tp@0 664 JJnumbofchars = ones(n1,1);
tp@0 665 if n1>9
tp@0 666 JJnumbofchars(10:n1) = JJnumbofchars(10:n1)+1;
tp@0 667 if n1 > 99
tp@0 668 JJnumbofchars(100:n1) = JJnumbofchars(100:n1)+1;
tp@0 669 if n1 > 999
tp@0 670 JJnumbofchars(1000:n1) = JJnumbofchars(1000:n1)+1;
tp@0 671 end
tp@0 672 end
tp@0 673 end
tp@0 674
tp@0 675 NSOU = int2str(nsources);
tp@0 676 NREC = int2str(nreceivers);
tp@0 677
tp@0 678 tsetup = etime(clock,t00);
tp@0 679
tp@0 680 % #########################################################################
tp@0 681 %
tp@0 682 % Find the paths
tp@0 683 %
tp@0 684 % #########################################################################
tp@0 685
tp@0 686 t00 = clock;
tp@0 687 tISEStree = zeros(nsources,1);
tp@0 688
tp@0 689 if calcpaths == 1
tp@0 690
tp@0 691 if SHOWTEXT >= 1
tp@0 692 disp(' ')
tp@0 693 disp('####################################################################')
tp@0 694 disp('#')
tp@0 695 disp('# Finding the paths')
tp@0 696 disp(' ')
tp@0 697 end
tp@0 698
tp@0 699 % Some variables that are needed for empty combinations
tp@0 700
tp@0 701 pathtypevec = [];
tp@0 702 reflpaths = [];
tp@0 703 pathlengthvec = [];
tp@0 704 specextradata = [];
tp@0 705 edgeextradata = [];
tp@0 706 mainlistguide = [];
tp@0 707 Sinsideplanenumber = [];
tp@0 708 Rinsideplanenumber = [];
tp@0 709 mainlistguide = [];
tp@0 710 mainlistguidepattern = [];
tp@0 711 directsoundrow = [];
tp@0 712 allspecrows = [];
tp@0 713 firstdiffrow = [];
tp@0 714 Varlist = [' pathtypevec reflpaths specextradata edgeextradata S R mainlistguide'];
tp@0 715 Varlist = [Varlist,' Sinsideplanenumber Rinsideplanenumber mainlistguide mainlistguidepattern directsoundrow allspecrows firstdiffrow'];
tp@0 716
tp@0 717 if doSbatches == 0
tp@0 718 eval(['load ',sdatafile])
tp@0 719 else
tp@0 720 latestsdatafile = reftoSbatchlist(1);
tp@0 721 [sdatafilepath,filestem,temp1] = fileparts(sdatafile);
tp@0 722 eval(['load ',[sdatafilepath,filesep,filestem],'_',int2str(latestsdatafile),'.mat'])
tp@0 723 end
tp@0 724 % % % % We create a souinsideplane matrix from the visplanesfroms matrix
tp@0 725 % % % % for use in EDB1makeirs
tp@0 726
tp@0 727 if doRbatches == 0
tp@0 728 eval(['load ',rdatafile])
tp@0 729 else
tp@0 730 latestrdatafile = reftoRbatchlist(1)
tp@0 731 [rdatafilepath,filestem,temp1] = fileparts(rdatafile)
tp@0 732 if ~isempty(rdatafilepath)
tp@0 733 txtstr = ['load ',[rdatafilepath,filesep,filestem],'_',int2str(latestrdatafile),'.mat']
tp@0 734 else
tp@0 735 txtstr = ['load ',[filestem],'_',int2str(latestrdatafile),'.mat']
tp@0 736 end
tp@0 737 eval(txtstr)
tp@0 738 end
tp@0 739 % % % % We create a recinsideplane matrix from the visplanesfromr matrix
tp@0 740 % % % % for use in EDB1makeirs
tp@0 741
tp@0 742 for isou = soulist_for_findpaths
tp@0 743
tp@0 744 t00_sou = clock;
tp@0 745 ISOU = int2str(isou);
tp@0 746 if SHOWTEXT >= 1
tp@0 747 disp(['Calculating for source ',ISOU,' (of ',NSOU,') '])
tp@0 748 end
tp@0 749 if doSbatches == 1
tp@0 750 if latestsdatafile ~= reftoSbatchlist(isou)
tp@0 751 latestsdatafile = reftoSbatchlist(isou);
tp@0 752 [sdatafilepath,filestem,temp1] = fileparts(sdatafile);
tp@0 753 eval(['load ',[sdatafilepath,filesep,filestem],'_',int2str(latestsdatafile),'.mat'])
tp@0 754 end
tp@0 755 Scolnumber = isou-Sbatchlist(reftoSbatchlist(isou),1)+1
tp@0 756 else
tp@0 757 Scolnumber = isou;
tp@0 758 end
tp@0 759 S = sources(Scolnumber,:);
tp@0 760
tp@0 761 % ###########################################################
tp@0 762 % #
tp@0 763 % # Calculate the ISES tree for each source
tp@0 764 % # IS = image sources, for specular reflections
tp@0 765 % # ES = edge sources, for edge diffraction
tp@0 766 % #
tp@0 767 % ###########################################################
tp@0 768
tp@0 769 if exist('ISEStreefile') ~= 1
tp@0 770 if SHOWTEXT >= 2
tp@0 771 disp([' Calculating ISES tree'])
tp@0 772 end
tp@0 773 if difforder >= 1
tp@0 774 usedISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms(:,isou),vispartedgesfroms(:,isou),nedgesubs);
tp@0 775 else
tp@0 776 usedISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms(:,isou),[],nedgesubs);
tp@0 777 end
tp@0 778
tp@0 779 else
tp@0 780 [ISEStreefilepath,filestem,temp1] = fileparts(ISEStreefile);
tp@0 781 usedISEStreefile = [[ISEStreefilepath,filesep,filestem],'_',ISOU,'_ISEStree.mat;'];
tp@0 782 if SHOWTEXT >= 2
tp@0 783 disp([' Using existing ISEStreefile: ',ISEStreefile,'_',ISOU,'_ISEStree.mat'])
tp@0 784 end
tp@0 785 end
tp@0 786 if SUPPRESSFILES == 0
tp@0 787 eval(['load ',usedISEStreefile])
tp@0 788 else
tp@0 789 POTENTIALISES = usedISEStreefile.POTENTIALISES;
tp@0 790 ORIGINSFROM = usedISEStreefile.ORIGINSFROM;
tp@0 791 ISCOORDS = usedISEStreefile.ISCOORDS;
tp@0 792 ISESVISIBILITY = usedISEStreefile.ISESVISIBILITY;
tp@0 793 IVNSPECMATRIX = usedISEStreefile. IVNSPECMATRIX;
tp@0 794 lengthNspecmatrix = usedISEStreefile.lengthNspecmatrix;
tp@0 795 IVNDIFFMATRIX = usedISEStreefile.IVNDIFFMATRIX;
tp@0 796 lengthNdiffmatrix = usedISEStreefile.lengthNdiffmatrix;
tp@0 797 singlediffcol = usedISEStreefile.singlediffcol;
tp@0 798 REFLORDER = usedISEStreefile.REFLORDER;
tp@0 799 startindicessinglediff = usedISEStreefile.startindicessinglediff;
tp@0 800 endindicessinglediff = usedISEStreefile.endindicessinglediff;
tp@0 801 ndecimaldivider = usedISEStreefile.ndecimaldivider;
tp@0 802 PointertoIRcombs = usedISEStreefile.PointertoIRcombs;
tp@0 803 IRoriginsfrom = usedISEStreefile.IRoriginsfrom;
tp@0 804 end
tp@0 805 tISEStree(isou) = etime(clock,t00_sou);
tp@0 806
tp@0 807 if dobackscatter ~= 1,
tp@0 808 reclist_in_forloop = reclist_for_findpaths;
tp@0 809 else
tp@0 810 reclist_in_forloop = isou;
tp@0 811 end
tp@0 812
tp@0 813 for irec = reclist_in_forloop
tp@0 814
tp@0 815 % ###########################################################
tp@0 816 % #
tp@0 817 % # Find the valid paths for each source-receiver combination
tp@0 818 % # _edpaths files are created
tp@0 819 % #
tp@0 820 % ###########################################################
tp@0 821
tp@0 822 IREC = int2str(irec);
tp@0 823 if SHOWTEXT >= 1
tp@0 824 if round(irec/idisp)*idisp==irec,
tp@0 825 disp(['...receiver no. ',IREC,' (of ',NREC,')'])
tp@0 826 end
tp@0 827 end
tp@0 828 if doRbatches == 1
tp@0 829 if latestrdatafile ~= reftoRbatchlist(irec)
tp@0 830 latestrdatafile = reftoRbatchlist(irec);
tp@0 831 [rdatafilepath,filestem,temp1] = fileparts(rdatafile);
tp@0 832 if ~isempty(rdatafilepath)
tp@0 833 txtstr = ['load ',[rdatafilepath,filesep,filestem],'_',int2str(latestrdatafile),'.mat'];
tp@0 834 else
tp@0 835 txtstr = ['load ',[filestem],'_',int2str(latestrdatafile),'.mat'];
tp@0 836 end
tp@0 837 eval(txtstr)
tp@0 838 end
tp@0 839 Rcolnumber = irec-Rbatchlist(reftoRbatchlist(irec),1)+1;
tp@0 840 else
tp@0 841 Rcolnumber = irec;
tp@0 842 end
tp@0 843 R = receivers(Rcolnumber,:);
tp@0 844
tp@0 845 if exist('edpathsfile') ~= 1
tp@0 846 desirededpathsfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0 847 if soutsidemodel(Scolnumber) == 0 & routsidemodel(Rcolnumber) == 0,
tp@0 848 if difforder > 0
tp@0 849 usededpathsfile = EDB1findpathsISESx(eddatafile,lengthNspecmatrix,lengthNdiffmatrix,singlediffcol,startindicessinglediff,...
tp@0 850 endindicessinglediff,ndecimaldivider,PointertoIRcombs,IRoriginsfrom,...
tp@0 851 S,R,Scolnumber,Rcolnumber,directsound,specorder,difforder,...
tp@0 852 nedgesubs,visplanesfroms(:,Scolnumber),visplanesfromr(:,Rcolnumber),...
tp@0 853 vispartedgesfroms(:,Scolnumber),vispartedgesfromr(:,Rcolnumber),desirededpathsfile);
tp@0 854 else
tp@0 855 usededpathsfile = EDB1findpathsISESx(eddatafile,lengthNspecmatrix,[],[],[],...
tp@0 856 [],[],[],[],S,R,Scolnumber,Rcolnumber,directsound,specorder,difforder,...
tp@0 857 nedgesubs,visplanesfroms(:,Scolnumber),visplanesfromr(:,Rcolnumber),[],[],desirededpathsfile);
tp@0 858 end
tp@0 859 else
tp@0 860 usededpathsfile = desirededpathsfile;
tp@0 861 eval(['save ',usededpathsfile,Varlist])
tp@0 862 end
tp@0 863 else
tp@0 864 if SHOWTEXT >= 2
tp@0 865 disp([' Using existing edpathsfile: ',edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'])
tp@0 866 desirededpathsfile = [edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0 867 end
tp@0 868 end
tp@0 869
tp@0 870
tp@0 871 end % of the for isou =
tp@0 872
tp@0 873 end % of the for irec =
tp@0 874 end
tp@0 875
tp@0 876 tfindpaths = etime(clock,t00);
tp@0 877
tp@0 878 % #########################################################################
tp@0 879 %
tp@0 880 % Edit paths
tp@0 881 %
tp@0 882 % #########################################################################
tp@0 883
tp@0 884 t00 = clock;
tp@0 885
tp@0 886 if ~isempty(symmetricedgepairs) & exist('edpathsfile') == 1
tp@0 887
tp@0 888 disp(['WARNING: You specified symmetric edge pairs, but path pruning is not'])
tp@0 889 disp([' allowed when you have specified an edpaths file'])
tp@0 890
tp@0 891 end
tp@0 892 if ~isempty(symmetricedgepairs) & exist('edpathsfile') ~= 1
tp@0 893
tp@0 894 if SHOWTEXT >= 1
tp@0 895 disp(' ')
tp@0 896 disp('####################################################################')
tp@0 897 disp('#')
tp@0 898 disp('# Pruning paths: removing redundant (symmetric) diffraction combinations')
tp@0 899 disp(' ')
tp@0 900 end
tp@0 901
tp@0 902 for isou = soulist_for_calcirs
tp@0 903 ISOU = int2str(isou);
tp@0 904 if SHOWTEXT >= 1
tp@0 905 disp(['Pruning paths for source ',ISOU,' (of ',NSOU,') '])
tp@0 906 end
tp@0 907
tp@0 908 if dobackscatter ~= 1,
tp@0 909 reclist_in_forloop = reclist_for_calcirs;
tp@0 910 else
tp@0 911 reclist_in_forloop = isou;
tp@0 912 end
tp@0 913
tp@0 914 for irec = reclist_in_forloop
tp@0 915 IREC = int2str(irec);
tp@0 916 if SHOWTEXT >= 1
tp@0 917 if round(irec/idisp)*idisp==irec,
tp@0 918 disp(['...receiver ',IREC,' (of ',NREC,')'])
tp@0 919 end
tp@0 920 end
tp@0 921
tp@0 922 edpathsfiletoedit = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0 923
tp@0 924 EDB1editpathlist(edpathsfiletoedit,[],symmetricedgepairs,[]);
tp@0 925
tp@0 926 end
tp@0 927 end
tp@0 928 end
tp@0 929
tp@0 930 tprunepaths = etime(clock,t00);
tp@0 931
tp@0 932 % #########################################################################
tp@0 933 %
tp@0 934 % Construct IRs
tp@0 935 %
tp@0 936 % #########################################################################
tp@0 937
tp@0 938 t00 = clock;
tp@0 939
tp@0 940 approxplanemidpoints = 0.5*(maxvals+minvals);
tp@0 941
tp@0 942 if calcirs == 1
tp@0 943
tp@0 944 if SHOWTEXT >= 1
tp@0 945 disp(' ')
tp@0 946 disp('####################################################################')
tp@0 947 disp('#')
tp@0 948 disp('# Constructing IRs')
tp@0 949 disp(' ')
tp@0 950 end
tp@0 951
tp@0 952 % If the parameter guiderowstouse wasn't specified in the setup file
tp@0 953 % we give it the default value.
tp@0 954
tp@0 955 if exist('guiderowstouse') ~= 1
tp@0 956 guiderowstouse = [];
tp@0 957 end
tp@0 958
tp@0 959 % Some variables that are needed for empty combinations
tp@0 960 irgeom = [];
tp@0 961 irtot = [];
tp@0 962 irdirect = [];
tp@0 963 irdiff = [];
tp@0 964 Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR elemsize nedgesubs'];
tp@0 965 Varlistaccum = [' irdirectaccum irgeomaccum irdiffaccum irtotaccum FSAMP Rstart CAIR elemsize nedgesubs'];
tp@0 966
tp@0 967 reclist_in_forloop = reclist_for_calcirs;
tp@0 968
tp@0 969 for irec = reclist_in_forloop
tp@0 970 IREC = int2str(irec);
tp@0 971 if SHOWTEXT >= 1
tp@0 972 disp(['...receiver ',IREC,' (of ',NREC,')'])
tp@0 973 end
tp@0 974
tp@0 975 irdirectaccum = [];
tp@0 976 irgeomaccum = [];
tp@0 977 irtotaccum = [];
tp@0 978 irdiffaccum = [];
tp@0 979
tp@0 980 if dobackscatter == 1,
tp@0 981 soulist_for_calcirs = irec;
tp@0 982 end
tp@0 983
tp@0 984 for isou = soulist_for_calcirs
tp@0 985 ISOU = int2str(isou);
tp@0 986 if SHOWTEXT >= 1
tp@0 987 if round(isou/idisp)*idisp==isou,
tp@0 988 disp(['Constructing IRs for source ',ISOU,' (of ',NSOU,') '])
tp@0 989 end
tp@0 990 end
tp@0 991
tp@0 992 if doaddsources == 1
tp@0 993 desiredirfile = [Filepath,Filestem,'_',IREC,'_ir.mat'];
tp@0 994 else
tp@0 995 desiredirfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_ir.mat'];
tp@0 996 end
tp@0 997
tp@0 998 if exist('edpathsfile') ~= 1
tp@0 999 usededpathsfile = [Filepath,Filestem,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0 1000 else
tp@0 1001 if SHOWTEXT >= 2
tp@0 1002 if round(irec/idisp)*idisp==irec,
tp@0 1003 disp([' Using existing edpathsfile: ',edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'])
tp@0 1004 end
tp@0 1005 end
tp@0 1006 usededpathsfile = [edpathsfile,'_',ISOU,'_',IREC,'_edpaths.mat'];
tp@0 1007 end
tp@0 1008
tp@0 1009 if difforder >= 2
tp@0 1010
tp@0 1011 edirfile = EDB1makeirs(usededpathsfile,specorder,...
tp@0 1012 Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0 1013 edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,...
tp@0 1014 reftoshortlistE,re1sho,re2sho,thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,directsound,saveindividualdiffirs);
tp@0 1015 else
tp@0 1016
tp@0 1017 edirfile = EDB1makeirs(usededpathsfile,specorder,...
tp@0 1018 Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0 1019 edgelengthvec,planeeqs,[],reflfactors,closwedangvec,planesatedge,elemsize,[],[],[],[],[],[],[],[],[],desiredirfile,guiderowstouse,directsound,saveindividualdiffirs);
tp@0 1020
tp@0 1021 end
tp@0 1022
tp@0 1023 if doaddsources == 1
tp@0 1024 eval(['load ',edirfile])
tp@0 1025
tp@0 1026 nirnew = length(irtot);
tp@0 1027 nirold = length(irtotaccum);
tp@0 1028 if nirnew > nirold
tp@0 1029 irtotaccum = [irtotaccum;zeros(nirnew-nirold,1)];
tp@0 1030 irdirectaccum = [irdirectaccum;zeros(nirnew-nirold,1)];
tp@0 1031 irgeomaccum = [irgeomaccum;zeros(nirnew-nirold,1)];
tp@0 1032 irdiffaccum = [irdiffaccum;zeros(nirnew-nirold,size(irdiff,2))];
tp@0 1033 end
tp@0 1034 irtotaccum(1:nirnew) = irtotaccum(1:nirnew) + irtot;
tp@0 1035 irdirectaccum(1:nirnew) = irdirectaccum(1:nirnew) + irdirect;
tp@0 1036 irgeomaccum(1:nirnew) = irgeomaccum(1:nirnew) + irgeom;
tp@0 1037 irdiffaccum(1:nirnew,:) = irdiffaccum(1:nirnew,:) + irdiff;
tp@0 1038 end
tp@0 1039
tp@0 1040 end % of the for isou = 1:nsou
tp@0 1041
tp@0 1042 if doaddsources == 1
tp@0 1043 eval(['save ',edirfile,Varlistaccum])
tp@0 1044 end
tp@0 1045
tp@0 1046 end % of the for irec = 1:nrec
tp@0 1047 end
tp@0 1048
tp@0 1049 %--------------------------------------------------------
tp@0 1050 % Save the timing data
tp@0 1051
tp@0 1052 tmakeirs = etime(clock,t00);
tp@0 1053
tp@0 1054 logfilename = [Filepath,Filestem,'_log.mat'];
tp@0 1055 Varlist = [' tcadgeo tedgeo tsgeo trgeo tsetup tISEStree tfindpaths tmakeirs tprunepaths nsources nreceivers'];
tp@0 1056 eval(['save ',logfilename,Varlist])