Mercurial > hg > human-echolocation

annotate private/EDB1makeirs.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 90220f7249fc
children
rev   line source
tp@0 1 function edirfile = EDB1makeirs(edpathsfile,specorder,Rstart,...
tp@0 2 EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0 3 edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,...
tp@0 4 thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,...
tp@0 5 saveindividualdiffirs)
tp@0 6 % EDB1makeirs - Constructs impulse responses from a list of paths in a file.
tp@0 7 %
tp@0 8 % Input parameters:
tp@0 9 % edpathsfile The name of the file containing the plane and edge data
tp@0 10 % specorder The maximum order of reflections that is calculated.
tp@0 11 % Rstart The reference distance that the time zero of the impulse
tp@0 12 % response refers to.
tp@0 13 % EDcalcmethod 'n' or 'v' or 'k', meaning:
tp@0 14 % 'n': the new method by Svensson et al
tp@0 15 % 'v': Vanderkooys method
tp@0 16 % 'k': the Kirchhoff diffraction approximation
tp@0 17 % edgestartcoords A matrix [nedges,3] of the startpoint coordinates
tp@0 18 % of
tp@0 19 % the edges.
tp@0 20 % edgeendcoords A matrix [nedges,3] of the endpoint coordinates of
tp@0 21 % the edges.
tp@0 22 % edgenvecs A matrix [nedges,3] of the normal vectors of
tp@0 23 % the reference plane of each edge.
tp@0 24 % edgelengthvec A list [nedge,1] of the lenghts of each edge.
tp@0 25 % planeeqs A matrix, [nplanes,4], of all the plane equations.
tp@0 26 % approxplanemidpoints A matrix, [nplanes,3], of all plane midpoint
tp@0 27 % coordinates. It may be an approximate location of
tp@0 28 % the midpoint.
tp@0 29 % reflfactors A list [nplanes,1] of the reflection factors of
tp@0 30 % each plane.
tp@0 31 % closwedangvec A list [nedge,1] of the close-wedge angle of each
tp@0 32 % edge.
tp@0 33 % planesatedge A matrix [nedge,2] of the two planes of each edge
tp@0 34 % with the reference plane first.
tp@0 35 % elemsize A list [1,difforder] with the relative edge element
tp@0 36 % sizes that will be used for each order. The first
tp@0 37 % value, for diffraction order 1, is redundant and
tp@0 38 % not used. For higher orders, a value of, e.g., 0.5
tp@0 39 % means that edges will be subdivided into elements
tp@0 40 % such that the wavelength at half the sampling
tp@0 41 % frequency is 0.5 times the element size. A lower
tp@0 42 % value gives faster but less accurate results.
tp@0 43 % reftoshortlistE A matrix, [nedges,nedges] with a pointer to the
tp@0 44 % short lists that contain edge-to-edge data.
tp@0 45 % re1sho A short list, [nshort,1] of the cylindrical radii
tp@0 46 % of each edge startpoint relative to all other
tp@0 47 % edges.
tp@0 48 % re2sho A short list, [nshort,1] of the cylindrical radii
tp@0 49 % of each edge endpoint relative to all other
tp@0 50 % edges.
tp@0 51 % thetae1sho A short list, [nshort,1] of the theta angle
tp@0 52 % of each edge startpoint relative to all other
tp@0 53 % edges.
tp@0 54 % thetae2sho A short list, [nshort,1] of the theta angle
tp@0 55 % of each edge endpoint relative to all other
tp@0 56 % edges.
tp@0 57 % ze1sho A short list, [nshort,1] of the z-value
tp@0 58 % of each edge startpoint relative to all other
tp@0 59 % edges.
tp@0 60 % ze2sho A short list, [nshort,1] of the z-value
tp@0 61 % of each edge endpoint relative to all other
tp@0 62 % edges.
tp@0 63 % edgeseespartialedge A matrix, [nedges,nedges], with edge-to-edge
tp@0 64 % visibility values. The values are 0 (invisible)
tp@0 65 % up to (2^(nedgesubs)-1)^2 (full visibility).
tp@0 66 % A pos. value indicates that the edge-to-edge path
tp@0 67 % runs along a plane; a neg. avlue indicates that it
tp@0 68 % does not run along a plane.
tp@0 69 % edgeplaneperptoplane1 A matrix, [nplanes,nedges], with 0 or 1. A
tp@0 70 % value 1 indicates that one of the edge's two
tp@0 71 % defining planes is perpendicular to another plane.
tp@0 72 % desiredirfile The file name that will be given to the
tp@0 73 % output file.
tp@0 74 % guiderowstouse (optional) A list of values, indicating which rows in the
tp@0 75 % mainlistguide and mainlistguidepattern that should
tp@0 76 % be used. This way it is possible to select only
tp@0 77 % diffraction, for instance. If this list is not
tp@0 78 % specified, all components will be used.
tp@0 79 % includedirectsound (optional) 0 or 1, indicating whether the direct
tp@0 80 % sound should be included or not. Default: 1
tp@0 81 % saveindividualdiffirs (optional) Vector of length 2 with values 0 or 1
tp@0 82 % Pos 1 indicating whether
tp@0 83 % 0 - all orders of diffraction IR:s will be summed
tp@0 84 % in a single vector 'irdiff', or
tp@0 85 % 1 - each order of diffraction IR:s will be placed
tp@0 86 % in individual columns in a matrix 'irdiff'.
tp@0 87 % Pos 2 indicating whether
tp@0 88 % 0 - only the total diffraction ir will be saved in the file.
tp@0 89 % 1 - all individual diffraction irs will be saved in a large
tp@0 90 % matrix alldiffirs.
tp@0 91 % Default: [0 0].
tp@0 92 % FSAMP, CAIR, SHOWTEXT Global parameters.
tp@0 93 %
tp@0 94 % Output parameters:
tp@0 95 % edirfile The filename used for the output file that contains
tp@0 96 % the impulse responses.
tp@0 97 %
tp@0 98 % Data in the output file:
tp@0 99 % irdirect The IR containing the direct sound, size [nirlength,1].
tp@0 100 % irgeom The IR containing the specular reflections, size [nirlength,1].
tp@0 101 % irdiff The IR containing the diffracted components, size [nirlength,1].
tp@0 102 % irtot The IR containing the sum of all components, size [nirlength,1].
tp@0 103 % alldiffirs (optional) Matrix with all individual first-order diffraction IRs
tp@0 104 % on individual lines.
tp@0 105 % alldiffirsextradata (optional) Vector with combination number (form reflpaths) that
tp@0 106 % matches alldiffirs.
tp@0 107 % FSAMP Rstart CAIR Same as the input parameters.
tp@0 108 %
tp@0 109 % Uses functions EDB1irfromslotlist EDB1calcdist EDB1coordtrans1 EDB1coordtrans2
tp@0 110 % EDB1wedge1st_int EDB1wedge2nd
tp@0 111 %
tp@0 112 % ----------------------------------------------------------------------------------------------
tp@0 113 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 114 %
tp@0 115 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 116 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 117 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 118 %
tp@0 119 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 120 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 121 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 122 %
tp@0 123 % You should have received a copy of the GNU General Public License along with the
tp@0 124 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 125 % ----------------------------------------------------------------------------------------------
tp@0 126 % Peter Svensson (svensson@iet.ntnu.no) 20061208
tp@0 127 %
tp@0 128 % edirfile = EDB1makeirs(edpathsfile,specorder,...
tp@0 129 % Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,...
tp@0 130 % edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,...
tp@0 131 % thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,...
tp@0 132 % saveindividualdiffirs);
tp@0 133
tp@0 134 global SHOWTEXT FSAMP CAIR BIGEDGESTEPMATRIX
tp@0 135
tp@0 136 global IRDIFFVEC
tp@0 137
tp@0 138 eval(['load ',edpathsfile])
tp@0 139
tp@0 140 %--------------------------------------------------------------------------
tp@0 141 % We look for the optional vector multfactors in the edpathsfile.
tp@0 142 % If it was there, then those values will be used to switch off, or boost
tp@0 143 % any reflection path.
tp@0 144 % Those values will only be used for diffraction paths.
tp@0 145
tp@0 146 if exist('multfactors') == 0
tp@0 147 multfactors = ones(size(reflpaths,1),1);
tp@0 148 else
tp@0 149 if size(multfactors,1) ~= size(reflpaths,1)
tp@0 150 error(['The edpaths file contained a vector multfactors which does not have the same length as reflpaths.'])
tp@0 151 end
tp@0 152 end
tp@0 153
tp@0 154 %--------------------------------------------------------------------------
tp@0 155
tp@0 156 edirfile = desiredirfile;
tp@0 157
tp@0 158 Sdirection = [1 0 0];
tp@0 159 maxnedgesorplanes = max(max(reflpaths));
tp@0 160 if ~isempty(maxnedgesorplanes)
tp@0 161 multfac = 10^ceil(log10(double(maxnedgesorplanes)+1));
tp@0 162 else
tp@0 163 multfac = 0;
tp@0 164 end
tp@0 165
tp@0 166 if nargin < 27
tp@0 167 saveindividualdiffirs = [0 0];
tp@0 168 end
tp@0 169 if nargin < 26
tp@0 170 includedirectsound = 1;
tp@0 171 end
tp@0 172 if nargin < 25
tp@0 173 guiderowstouse = [];
tp@0 174 end
tp@0 175 if isempty(guiderowstouse)
tp@0 176 usesubset = 0;
tp@0 177 if SHOWTEXT > 2
tp@0 178 disp('Using all components')
tp@0 179 end
tp@0 180 else
tp@0 181 usesubset = 1;
tp@0 182 if SHOWTEXT > 2
tp@0 183 disp('Using only some components')
tp@0 184 end
tp@0 185 end
tp@0 186
tp@0 187 nyvec = pi./(2*pi - closwedangvec);
tp@0 188 onesvec1 = ones(1,3);
tp@0 189
tp@0 190 souspecboost = 1;
tp@0 191 if ~isempty(Sinsideplanenumber)
tp@0 192 if reflfactors(Sinsideplanenumber(1)) ~= 0
tp@0 193 souspecboost = 2;
tp@0 194 end
tp@0 195 end
tp@0 196
tp@0 197 recspecboost = 1;
tp@0 198 if ~isempty(Rinsideplanenumber)
tp@0 199 if reflfactors(Rinsideplanenumber(1)) ~= 0
tp@0 200 recspecboost = 2;
tp@0 201 end
tp@0 202 end
tp@0 203
tp@0 204 if isempty(re1sho)
tp@0 205 userwantsdiff2 = 0;
tp@0 206 else
tp@0 207 userwantsdiff2 = 1;
tp@0 208 end
tp@0 209
tp@0 210 %-------------------------------------------------------
tp@0 211
tp@0 212 if exist('mainlistguide') ~= 1
tp@0 213 mainlistguide = [];
tp@0 214 else
tp@0 215 mainlistguide = double(mainlistguide);
tp@0 216 end
tp@0 217
tp@0 218 [ncomponents,ncols] = size(reflpaths);
tp@0 219 [nrowsinlist,slask] = size(mainlistguide);
tp@0 220 if ncols > 1
tp@0 221 difforderinlist = sum(mainlistguidepattern.'=='d').';
tp@0 222 else
tp@0 223 difforderinlist = (mainlistguidepattern=='d');
tp@0 224 end
tp@0 225
tp@0 226 lastNdiffrow = zeros(1,specorder);
tp@0 227 for ii = 1:specorder
tp@0 228 iv = find(difforderinlist <= ii );
tp@0 229 if ~isempty(iv)
tp@0 230 lastNdiffrow(ii) = iv(end);
tp@0 231 end
tp@0 232 end
tp@0 233
tp@0 234 nrefltogetherwdiff = sum(mainlistguidepattern.'=='d').'.*( sum(mainlistguidepattern.'=='s').');
tp@0 235
tp@0 236 if SHOWTEXT >= 3
tp@0 237 disp([' Constructing IR. ',int2str(ncomponents),' components.'])
tp@0 238 end
tp@0 239
tp@0 240 irdirect = [0 0].';
tp@0 241 irgeom = [0 0].';
tp@0 242 irdiff = [0 0].';
tp@0 243 irtot = [0 0].';
tp@0 244
tp@0 245 Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod'];
tp@0 246
tp@0 247 if ncomponents == 0
tp@0 248 eval(['save ',edirfile,Varlist])
tp@0 249 return
tp@0 250 end
tp@0 251
tp@0 252 %##############################################################################
tp@0 253 %##############################################################################
tp@0 254 %##############################################################################
tp@0 255 %##############################################################################
tp@0 256 %
tp@0 257 % Diffraction once, possibly with pre- and post-specular reflections
tp@0 258 %
tp@0 259 %##############################################################################
tp@0 260
tp@0 261 directsoundonboundary = 0;
tp@0 262 specreflonboundary = zeros(size(planeeqs,1),1);
tp@0 263
tp@0 264 if firstdiffrow ~= 0
tp@0 265
tp@0 266 % First we remove the 'd' for all the rows in the mainlistguidepattern that
tp@0 267 % we do not want to use.
tp@0 268
tp@0 269 if usesubset == 1
tp@0 270 singdiffrows = find(difforderinlist==1);
tp@0 271 rowstoclear = setdiff(singdiffrows,guiderowstouse);
tp@0 272 if ~isempty(rowstoclear)
tp@0 273 mainlistguidepattern(rowstoclear,:) = mainlistguidepattern(rowstoclear,:)*0;
tp@0 274 end
tp@0 275 end
tp@0 276
tp@0 277 % Here we can use the single-diffraction calculation for all
tp@0 278 % combinations. The source should either be the original source or
tp@0 279 % an image source. The receiver should either be the original
tp@0 280 % receiver or an image receiver.
tp@0 281
tp@0 282 % First we create lists that specify whether there are specular
tp@0 283 % reflections before the edge diffraction, and after the edge
tp@0 284 % diffraction. Also, a list which gives the column where the edge
tp@0 285 % diffraction is present can be used to efficiently extract the
tp@0 286 % edge number.
tp@0 287 % NB! The lists prespeclist, postspeclist and singdiffcol all have
tp@0 288 % the (short) length of the mainlistguide.
tp@0 289 % The data in these lists will then be used to create the long
tp@0 290 % lists prespecpresent and postspecpresent which simply are 0 or 1.
tp@0 291 % A matrix called edgemask will have the same number of columns as
tp@0 292 % the reflpaths matrix. For each row, there will be a 1 in the
tp@0 293 % column where the edge diffraction is so we can efficiently find
tp@0 294 % the edge number.
tp@0 295
tp@0 296 multmatrix = [1:specorder];
tp@0 297 if size(mainlistguidepattern,2) > specorder
tp@0 298 multmatrix = [1:size(mainlistguidepattern,2)];
tp@0 299 end
tp@0 300 multmatrix = multmatrix( ones(nrowsinlist,1),: );
tp@0 301
tp@0 302 if ncols > 1
tp@0 303 singdiffcol = sum( (multmatrix.*(mainlistguidepattern=='d')).' ).';
tp@0 304 singdiffcol = singdiffcol.*(sum(mainlistguidepattern.'=='d').'<=1);
tp@0 305 nrefl = sum( (mainlistguidepattern=='s' | mainlistguidepattern=='d').' ).';
tp@0 306 else
tp@0 307 singdiffcol = mainlistguidepattern=='d';
tp@0 308 nrefl = singdiffcol | mainlistguidepattern=='s';
tp@0 309 end
tp@0 310
tp@0 311 prespeclist = singdiffcol-1;
tp@0 312 prespeclist = prespeclist.*(prespeclist>0);
tp@0 313 postspeclist = (nrefl-singdiffcol).*(singdiffcol~=0);
tp@0 314
tp@0 315 prespecpresent = zeros(ncomponents,1);
tp@0 316 postspecpresent = zeros(ncomponents,1);
tp@0 317 edgemask = zeros(ncomponents,specorder);
tp@0 318
tp@0 319 diffrowsthatareOK = [firstdiffrow:lastNdiffrow(1)];
tp@0 320 if usesubset == 1
tp@0 321 diffrowsthatareOK = intersect(diffrowsthatareOK,guiderowstouse);
tp@0 322 end
tp@0 323 for ii = diffrowsthatareOK, %firstdiffrow:lastNdiffrow(1)
tp@0 324 iv = [(mainlistguide(ii,2)):(mainlistguide(ii,3))];
tp@0 325 onesvec2 = ones(length(iv),1);
tp@0 326 edgemask(iv,singdiffcol(ii)) = onesvec2;
tp@0 327 prespecpresent(iv) = (prespeclist(ii)>0)*onesvec2;
tp@0 328 postspecpresent(iv) = (postspeclist(ii)>0)*onesvec2;
tp@0 329 end
tp@0 330 prespecpresent = prespecpresent(:,onesvec1);
tp@0 331 postspecpresent = postspecpresent(:,onesvec1);
tp@0 332
tp@0 333 if ncols > 1
tp@0 334 edgenumb = sum( (edgemask.*double(reflpaths(:,1:specorder))).').';
tp@0 335 else
tp@0 336 edgenumb = edgemask.*double(reflpaths);
tp@0 337 end
tp@0 338
tp@0 339 nedgeirs = size(edgenumb,1);
tp@0 340 nnonzeroedgeirs = sum(edgenumb(:,1)>0);
tp@0 341
tp@0 342 ircounter = 0;
tp@0 343 edgelist = unique(edgenumb);
tp@0 344 for ii = 1: length(edgelist)
tp@0 345 edge = edgelist(ii);
tp@0 346 if edge~= 0
tp@0 347 iv = find(edgenumb==edge);
tp@0 348 ncombs = length(iv);
tp@0 349
tp@0 350 if ncombs > 0 & any(multfactors(iv))
tp@0 351 if SHOWTEXT >= 4
tp@0 352 disp([' Edge ',int2str(edge),': ',int2str(ncombs),' combinations'])
tp@0 353 end
tp@0 354
tp@0 355 onesvec3 = ones(ncombs,1);
tp@0 356
tp@0 357 IS = full(specextradata(iv,1:3));
tp@0 358 IR = full(specextradata(iv,4:6));
tp@0 359
tp@0 360 edgestart = full(edgeextradata(iv,1));
tp@0 361 edgeend = full(edgeextradata(iv,2));
tp@0 362
tp@0 363 % Calculate the cyl coordinates of all IS
tp@0 364 %
tp@0 365 % Calculate the cyl coord of all IR
tp@0 366
tp@0 367 edgecoords = [edgestartcoords(edge,:);edgeendcoords(edge,:)];
tp@0 368 [rs,thetas,zs,rr,thetar,zr] = EDB1coordtrans2(IS,IR,edgecoords,edgenvecs(edge,:));
tp@0 369
tp@0 370 bc = reflfactors(planesatedge(edge,:)).';
tp@0 371 if prod(double(bc==[1 1])) ~= 1
tp@0 372 disp(' Non-rigid wedge')
tp@0 373 end
tp@0 374
tp@0 375 for jj = 1:ncombs
tp@0 376 if multfactors(iv(jj)) > 0
tp@0 377 if EDcalcmethod(1) == 'n'
tp@0 378
tp@0 379 [irnew,slask,singularterm] = EDB1wedge1st_int(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),...
tp@0 380 edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc);
tp@0 381
tp@0 382 if any(singularterm)
tp@0 383 if singularterm(2) | singularterm(3)
tp@0 384 directsoundonboundary = 1;
tp@0 385 elseif singularterm(1)
tp@0 386 if specorder == 1
tp@0 387 specreflonboundary(planesatedge(edge,2)) = 1;
tp@0 388 else
tp@0 389 disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!');
tp@0 390 end
tp@0 391 elseif singularterm(4)
tp@0 392 if specorder == 1
tp@0 393 specreflonboundary(planesatedge(edge,1)) = 1;
tp@0 394 else
tp@0 395 disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!');
tp@0 396 end
tp@0 397 end
tp@0 398 end
tp@0 399
tp@0 400 else % ... if EDcalcmethod(1) == 'n'
tp@0 401
tp@0 402 [irnew,slask,singularterm] = EDB1wedge1stcombo(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),...
tp@0 403 edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc);
tp@0 404
tp@0 405 end % ... if EDcalcmethod(1) == 'n'
tp@0 406
tp@0 407 % Decide whether the IR should be boosted or not
tp@0 408 %
tp@0 409 % The factors souspecboost and recspecboost have
tp@0 410 % values other than 1 if the source or receiver is directly
tp@0 411 % at a plane, for the boosting of the direct sound, and
tp@0 412 % specular reflections.
tp@0 413 %
tp@0 414 % This boost factor should be used also for edge
tp@0 415 % diffraction if:
tp@0 416 % 1. There are some specular reflections between the
tp@0 417 % source and the edge
tp@0 418 % or
tp@0 419 % 2. There are no specular reflections between the
tp@0 420 % source and the edge, but the source/receiver is at
tp@0 421 % a plane which does not belong to the edge.
tp@0 422
tp@0 423 if souspecboost ~= 1 | recspecboost ~= 1
tp@0 424
tp@0 425 boostfactor = 1;
tp@0 426 if prespecpresent(iv(jj)) == 1
tp@0 427 boostfactor = souspecboost;
tp@0 428 else
tp@0 429 if ~isempty(Sinsideplanenumber)
tp@0 430 if Sinsideplanenumber(1)~=planesatedge(edge,1) & Sinsideplanenumber(1)~=planesatedge(edge,2),
tp@0 431 boostfactor = souspecboost;
tp@0 432 end
tp@0 433 end
tp@0 434 end
tp@0 435 if postspecpresent(iv(jj)) == 1
tp@0 436 boostfactor = boostfactor*recspecboost;
tp@0 437 else
tp@0 438 if ~isempty(Rinsideplanenumber)
tp@0 439 if Rinsideplanenumber(1)~=planesatedge(edge,1) & Rinsideplanenumber(1)~=planesatedge(edge,2),
tp@0 440 boostfactor = boostfactor*recspecboost;
tp@0 441 end
tp@0 442 end
tp@0 443 end
tp@0 444 if boostfactor ~= 1
tp@0 445 irnew = irnew*boostfactor;
tp@0 446 end
tp@0 447
tp@0 448 end % .... if souspecboost ~= 1 | recspecboost ~= 1
tp@0 449
tp@0 450 ndiff = length(irdiff);
tp@0 451 nnew = length(irnew);
tp@0 452
tp@0 453 ircounter = ircounter + 1;
tp@0 454
tp@0 455 if nnew > ndiff
tp@0 456 irdiff = [irdiff;zeros(nnew-ndiff,1)];
tp@0 457 end
tp@0 458 irdiff(1:nnew) = irdiff(1:nnew) + irnew*double(multfactors(iv(jj)));
tp@0 459
tp@0 460 if saveindividualdiffirs(2) == 1
tp@0 461 if exist('alldiffirs','var') == 0
tp@0 462 alldiffirs = sparse(zeros(nnonzeroedgeirs,nnew));
tp@0 463 alldiffirs(1,:) = irnew*double(multfactors(iv(jj))).';
tp@0 464 alldiffirsextradata = zeros(nnonzeroedgeirs,1);
tp@0 465 alldiffirsextradata(1) = iv(jj);
tp@0 466 else
tp@0 467 if nnew > ndiff
tp@0 468 alldiffirs = [alldiffirs zeros(nnonzeroedgeirs,nnew-ndiff)];
tp@0 469 end
tp@0 470 alldiffirs(ircounter,1:nnew) = irnew*double(multfactors(iv(jj))).';
tp@0 471 alldiffirsextradata(ircounter) = iv(jj);
tp@0 472 end
tp@0 473
tp@0 474 end
tp@0 475
tp@0 476
tp@0 477 if SHOWTEXT >= 7
tp@0 478 figure(1)
tp@0 479 plot(irnew)
tp@0 480 figure(2)
tp@0 481 plot(irdiff)
tp@0 482 save ~/Documents/Temp/irdiff2.mat
tp@0 483 pause
tp@0 484 end
tp@0 485
tp@0 486 end % ... if multfactors(iv) > 0
tp@0 487
tp@0 488 end % ..... for jj = 1:ncombs
tp@0 489
tp@0 490 end %.... if ncombs > 0 & any(multfactors(iv))
tp@0 491
tp@0 492 end % .... if edge~= 0
tp@0 493
tp@0 494 end % .... for ii = 1: length(edgelist)
tp@0 495
tp@0 496 end % ...if firstdiffrow ~= 0
tp@0 497
tp@0 498 nspecreflonboundary = sum(specreflonboundary);
tp@0 499
tp@0 500 %##############################################################################
tp@0 501 %##############################################################################
tp@0 502 %##############################################################################
tp@0 503 %##############################################################################
tp@0 504 %
tp@0 505 % The direct sound
tp@0 506 %
tp@0 507 %##############################################################################
tp@0 508
tp@0 509 if usesubset == 1 & directsoundrow == 1
tp@0 510 if any(guiderowstouse == 1) == 0
tp@0 511 directsoundrow = 0;
tp@0 512 end
tp@0 513 end
tp@0 514
tp@0 515 if directsoundrow == 1 & includedirectsound == 1
tp@0 516 dist = norm(R-S);
tp@0 517 slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1;
tp@0 518 amp = souspecboost*recspecboost./dist;
tp@0 519 if directsoundonboundary
tp@0 520 amp = amp/2;
tp@0 521 if SHOWTEXT >= 4
tp@0 522 disp('HALVING DIRECT SOUND')
tp@0 523 end
tp@0 524 end
tp@0 525 slotnumber = floor(slotnumberfrac);
tp@0 526 if any(slotnumber<1)
tp@0 527 error('ERROR: Rstart is set to too low a value!')
tp@0 528 end
tp@0 529 slotnumberfrac = slotnumberfrac - slotnumber;
tp@0 530 irdirect = zeros( slotnumber+2,1 );
tp@0 531
tp@0 532 irdirect(slotnumber) = amp.*(1-slotnumberfrac) ;
tp@0 533 irdirect(slotnumber+1) = amp.*slotnumberfrac;
tp@0 534
tp@0 535 if ncomponents == 1
tp@0 536 nnew = length(irdirect);
tp@0 537 if nnew > 2
tp@0 538 irdiff = [irdiff;zeros(nnew-2,1)];
tp@0 539 irtot = [irtot;zeros(nnew-2,1)];
tp@0 540 irgeom = [irgeom;zeros(nnew-2,1)];
tp@0 541 end
tp@0 542 irtot(1:nnew) = irtot(1:nnew) + irdirect;
tp@0 543
tp@0 544 irtot = sparse(irtot);
tp@0 545 irgeom = sparse(irgeom);
tp@0 546 irdirect = sparse(irdirect);
tp@0 547 irdiff = sparse(irdiff);
tp@0 548
tp@0 549 eval(['save ',edirfile,Varlist])
tp@0 550 return
tp@0 551 end
tp@0 552 else
tp@0 553 if directsoundonboundary == 1
tp@0 554
tp@0 555 dist = norm(R-S);
tp@0 556 slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1;
tp@0 557 amp = souspecboost*recspecboost./dist;
tp@0 558 if directsoundonboundary
tp@0 559 amp = amp/2;
tp@0 560 if SHOWTEXT >= 4
tp@0 561 disp('HALVING DIRECT SOUND')
tp@0 562 end
tp@0 563 end
tp@0 564 slotnumber = floor(slotnumberfrac);
tp@0 565 if any(slotnumber<1)
tp@0 566 error('ERROR: Rstart is set to too low a value!')
tp@0 567 end
tp@0 568 slotnumberfrac = slotnumberfrac - slotnumber;
tp@0 569 irdirect = zeros( slotnumber+2,1 );
tp@0 570
tp@0 571 irdirect(slotnumber) = amp.*(1-slotnumberfrac) ;
tp@0 572 irdirect(slotnumber+1) = amp.*slotnumberfrac;
tp@0 573
tp@0 574 end
tp@0 575 if SHOWTEXT >= 4
tp@0 576 disp([' No direct sound'])
tp@0 577 end
tp@0 578 end
tp@0 579
tp@0 580 %##############################################################################
tp@0 581 %##############################################################################
tp@0 582 %##############################################################################
tp@0 583 %##############################################################################
tp@0 584 %
tp@0 585 % All-specular s, ss, sss, etc
tp@0 586 %
tp@0 587 %##############################################################################
tp@0 588
tp@0 589 if allspecrows(1) ~= 0
tp@0 590 if usesubset == 1
tp@0 591 specrowsthatareOK = intersect(allspecrows,guiderowstouse);
tp@0 592 else
tp@0 593 specrowsthatareOK = [allspecrows(1):allspecrows(2)];
tp@0 594 end
tp@0 595 if ~isempty(specrowsthatareOK)
tp@0 596 if usesubset == 1
tp@0 597 temp = mainlistguide(specrowsthatareOK,2:3);
tp@0 598 ivspec = [double(mainlistguide(specrowsthatareOK(1),2)):double(mainlistguide(specrowsthatareOK(1),3))];
tp@0 599 for ii = 2:size(temp,1);
tp@0 600 ivspec = [ivspec [double(mainlistguide(specrowsthatareOK(ii),2)):double(mainlistguide(specrowsthatareOK(ii),3))]];
tp@0 601 end
tp@0 602 else
tp@0 603 ivspec = [mainlistguide(allspecrows(1),2):mainlistguide(allspecrows(2),3)];
tp@0 604 end
tp@0 605
tp@0 606 if nspecreflonboundary > 0
tp@0 607 listofspecreflonboundary = find(specreflonboundary);
tp@0 608 specrefltoadd = setdiff(listofspecreflonboundary,ivspec)
tp@0 609 if ~isempty(specrefltoadd)
tp@0 610 error(['ERROR: We need to add some code here, to reinsert pruned spec. refl.'])
tp@0 611 end
tp@0 612 end
tp@0 613
tp@0 614 if SHOWTEXT >= 4
tp@0 615 disp([' ',int2str(length(ivspec)),' specular reflections'])
tp@0 616 end
tp@0 617
tp@0 618 dist = EDB1calcdist(full(specextradata(ivspec,1:3)),R);
tp@0 619
tp@0 620 if nspecreflonboundary > 0
tp@0 621 specamp = ones(size(souspecboost));
tp@0 622 amplitudeshouldbehalf = ismember(reflpaths(ivspec,1),listofspecreflonboundary);
tp@0 623 specamp = specamp - amplitudeshouldbehalf*0.5;
tp@0 624 irnew = specamp*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);
tp@0 625 else
tp@0 626 irnew = souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);
tp@0 627 end
tp@0 628 ngeom = length(irgeom);
tp@0 629 nnew = length(irnew);
tp@0 630 if nnew > ngeom
tp@0 631 irgeom = [irgeom;zeros(nnew-ngeom,1)];
tp@0 632 end
tp@0 633 irgeom(1:nnew) = irgeom(1:nnew) + irnew;
tp@0 634
tp@0 635 else
tp@0 636 if nspecreflonboundary > 0
tp@0 637 listofspecreflonboundary = find(specreflonboundary);
tp@0 638 error(['ERROR: We need to add some code here, to reinsert pruned spec. refl., pos. 2'])
tp@0 639 end
tp@0 640 end
tp@0 641 else
tp@0 642 if nspecreflonboundary > 0
tp@0 643 listofspecreflonboundary = find(specreflonboundary);
tp@0 644
tp@0 645 [xis] = EDB1findis(S,listofspecreflonboundary,planeeqs,1,[1 1 1]);
tp@0 646
tp@0 647 disp(['WARNING: We have some new code here, to reinsert pruned spec. refl., pos. 3'])
tp@0 648 dist = EDB1calcdist(xis,R);
tp@0 649 irnew = 0.5*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist);
tp@0 650 ngeom = length(irgeom);
tp@0 651 nnew = length(irnew);
tp@0 652 if nnew > ngeom
tp@0 653 irgeom = [irgeom;zeros(nnew-ngeom,1)];
tp@0 654 end
tp@0 655 irgeom(1:nnew) = irgeom(1:nnew) + irnew;
tp@0 656
tp@0 657 end
tp@0 658 end
tp@0 659
tp@0 660 %##############################################################################
tp@0 661 %##############################################################################
tp@0 662 %##############################################################################
tp@0 663 %##############################################################################
tp@0 664 %
tp@0 665 % Multiple diffraction, possibly with pre- and post-specular reflections
tp@0 666 %
tp@0 667 %##############################################################################
tp@0 668
tp@0 669 if userwantsdiff2 == 1
tp@0 670
tp@0 671 JJ = setstr(32*ones(specorder,1));
tp@0 672 for jj=1:specorder
tp@0 673 jjstr = int2str(jj);
tp@0 674 JJ(jj,1:length(jjstr)) = jjstr;
tp@0 675 end
tp@0 676 [n1,n2] = size(JJ);
tp@0 677
tp@0 678
tp@0 679 for Ndifforder = 2:specorder
tp@0 680 if SHOWTEXT >= 3
tp@0 681 disp([' Diffraction order ',int2str(Ndifforder)])
tp@0 682 end
tp@0 683
tp@0 684 if any(difforderinlist==Ndifforder) & elemsize(Ndifforder) > 0
tp@0 685
tp@0 686 % Calculate some general parameters that are shared for all
tp@0 687 % N-diffraction calculations
tp@0 688
tp@0 689 divmin = CAIR/(FSAMP*elemsize(Ndifforder));
tp@0 690 ndivvec = ceil(abs( edgelengthvec.' )/divmin);
tp@0 691 dzvec = (edgelengthvec.')./ndivvec;
tp@0 692
tp@0 693 ncylrows = 4*(Ndifforder-1);
tp@0 694
tp@0 695 % Here we can use the double-diffraction calculation for all
tp@0 696 % combinations. The source should either be the original source or
tp@0 697 % an image source. The receiver should either be the original
tp@0 698 % receiver or an image receiver.
tp@0 699
tp@0 700 noffset = lastNdiffrow(Ndifforder-1);
tp@0 701 ivNdiff = [noffset+1:lastNdiffrow(Ndifforder)];
tp@0 702 ndiffcombs = length(ivNdiff);
tp@0 703 zerosvec1 = zeros(ndiffcombs,1);
tp@0 704
tp@0 705 ndoubrows = length(ivNdiff);
tp@0 706 previousrow = lastNdiffrow(Ndifforder-1);
tp@0 707 if previousrow > 0
tp@0 708 noffsetlonglist = mainlistguide(previousrow,3);
tp@0 709 else
tp@0 710 noffsetlonglist = 0;
tp@0 711 end
tp@0 712 nremaining = ncomponents - (mainlistguide(lastNdiffrow(Ndifforder),3));
tp@0 713 nlonglist = ncomponents-noffsetlonglist-nremaining;
tp@0 714 ivlonglist = [mainlistguide(ivNdiff(1),2):mainlistguide(ivNdiff(end),3)].';
tp@0 715
tp@0 716 diffcols = zerosvec1(:,ones(1,Ndifforder));
tp@0 717 for ii = ivNdiff(1):ivNdiff(end)
tp@0 718 diffcols(ii-ivNdiff(1)+1,:) = find(mainlistguidepattern(ii,:)=='d');
tp@0 719 end
tp@0 720
tp@0 721 nreflorder = sum( (mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='d'|mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='s').').';
tp@0 722 nprespecs = diffcols(:,1)-1;
tp@0 723 nmidspecs = diffcols(:,2)-diffcols(:,1)-1;
tp@0 724 npostspecs = nreflorder-diffcols(:,Ndifforder);
tp@0 725
tp@0 726 if ndiffcombs > 1
tp@0 727 ivreftolonglist = ivlonglist(:,ones(1,ndiffcombs));
tp@0 728 comppattern = mainlistguide(ivNdiff,2).';
tp@0 729 comppattern = comppattern(ones(nlonglist,1),:);
tp@0 730
tp@0 731 rowinpatternlist = sum( (ivreftolonglist>=comppattern).' ).';
tp@0 732 else
tp@0 733 rowinpatternlist = ones(size(ivlonglist));
tp@0 734 end
tp@0 735
tp@0 736 % Construct a long matrix with the edge numbers extracted from the
tp@0 737 % reflpaths matrix.
tp@0 738
tp@0 739 longnmidspec = zeros(nlonglist,1);
tp@0 740 iv = find(nmidspecs>0);
tp@0 741 for ii = 1:length(iv)
tp@0 742 ivreftolonglist = [mainlistguide(ivNdiff(iv(ii)),2):mainlistguide(ivNdiff(iv(ii)),3)] - noffsetlonglist;
tp@0 743 longnmidspec(ivreftolonglist) = nmidspecs(iv(ii));
tp@0 744 end
tp@0 745
tp@0 746 %------------------------------------------------------------------
tp@0 747 % edgepattern will be a matrix which, for each reflection combination, contains
tp@0 748 % the 2,3,4,... edge numbers that are involved in each path
tp@0 749 % regardless of there are specular reflections before, in the
tp@0 750 % middle, or after.
tp@0 751
tp@0 752 edgepattern = zeros(nlonglist,Ndifforder);
tp@0 753 countvec = [1:nlonglist].';
tp@0 754 reflpathscut = reflpaths(ivlonglist,:);
tp@0 755 for ii = 1:Ndifforder
tp@0 756 ivrefcol = countvec + (diffcols(rowinpatternlist,ii)-1)*nlonglist;
tp@0 757 edgepattern(:,ii) = reflpathscut(ivrefcol);
tp@0 758 end
tp@0 759
tp@0 760 %------------------------------------------------------------------
tp@0 761 % midspecpattern will be a matrix which, for each reflection combination, contains
tp@0 762 % the 1,2,3,4,... specular reflection numbers (i.e., plane numbers) that are involved
tp@0 763 % in each path between diffraction 1 and diffraction 2.
tp@0 764 % First, the reflection component immediately after diffraction 1
tp@0 765 % is selected for every reflection. Then there is a masking which
tp@0 766 % zeroes the midspecpattern value for all the combinations that
tp@0 767 % don't have any midspec.
tp@0 768 % NB! For combinations like ssssdd, the first step will point to a
tp@0 769 % column which is outside the matrix. This is fixed by maximizing
tp@0 770 % the column number.
tp@0 771
tp@0 772 midspecpattern = zeros(nlonglist,max([specorder-Ndifforder 1]));
tp@0 773 maxpointvalue = nlonglist*max([specorder-Ndifforder 1]);
tp@0 774 for ii = 1:specorder-Ndifforder
tp@0 775 ivrefcol = countvec + (diffcols(rowinpatternlist,1)+(ii-1))*nlonglist;
tp@0 776 ivrefcol = mod(ivrefcol-1,maxpointvalue)+1;
tp@0 777 midspecpattern(:,ii) = double(reflpathscut(ivrefcol)).*(longnmidspec>=ii);
tp@0 778 end
tp@0 779 diff1col = diffcols(rowinpatternlist,1);
tp@0 780
tp@0 781 % Many combinations include the same edge-pair/N-let, so we extract the
tp@0 782 % unique edge pairs/N-lets, and go through them in the ii-loop
tp@0 783
tp@0 784 [edgeshortlist,ivec,jvec] = unique(edgepattern,'rows');
tp@0 785 lastndivcomb = zeros(1,Ndifforder);
tp@0 786
tp@0 787 for ii = 1: length(ivec)
tp@0 788
tp@0 789 if SHOWTEXT >= 3
tp@0 790 if round(ii/ceil(length(ivec)/10))*ceil(length(ivec)/10) == ii
tp@0 791 disp([' Combination no. ',int2str(ii),' of ',int2str(length(ivec))])
tp@0 792 end
tp@0 793 end
tp@0 794
tp@0 795 iv = find(jvec==ii);
tp@0 796 ncombs = length(iv);
tp@0 797 onesvec3 = ones(ncombs,1);
tp@0 798
tp@0 799 if SHOWTEXT >= 4
tp@0 800 numvec = int2str(edgeshortlist(ii,1));
tp@0 801 for ll = 2:Ndifforder
tp@0 802 numvec = [numvec,' ', int2str(edgeshortlist(ii,ll))];
tp@0 803 end
tp@0 804 disp([' Edges ',numvec])
tp@0 805
tp@0 806 end
tp@0 807
tp@0 808 if Ndifforder >= 2 & any(multfactors(ivlonglist(iv)))
tp@0 809
tp@0 810 newndivvec = ndivvec(edgeshortlist(ii,:));
tp@0 811
tp@0 812 if any(lastndivcomb~=newndivvec)
tp@0 813 if Ndifforder == 2
tp@0 814 ivmatrix = EDB1creindexmatrix(newndivvec);
tp@0 815 else
tp@0 816 ivmatrix = EDB1creindexmatrix(newndivvec(2:end));
tp@0 817 end
tp@0 818 [nedgeelcombs,slask] = size(ivmatrix);
tp@0 819 if Ndifforder == 2
tp@0 820 BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),:);
tp@0 821 else
tp@0 822 BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),2:end);
tp@0 823 end
tp@0 824 clear ivmatrix
tp@0 825 lastndivcomb = newndivvec;
tp@0 826
tp@0 827 end
tp@0 828 end
tp@0 829
tp@0 830 pathalongplane = (edgeseespartialedge(edgeshortlist(ii,2),edgeshortlist(ii,1))>0);
tp@0 831 for jj = 3:Ndifforder
tp@0 832 pathalongplane = [pathalongplane,(edgeseespartialedge(edgeshortlist(ii,jj),edgeshortlist(ii,jj-1))>0)];
tp@0 833 end
tp@0 834
tp@0 835 if any(multfactors(ivlonglist(iv)))
tp@0 836
tp@0 837 IS = full(specextradata(ivlonglist(iv),1:3));
tp@0 838 IR = full(specextradata(ivlonglist(iv),4:6));
tp@0 839
tp@0 840 firstedgestart = full(edgeextradata(ivlonglist(iv),1));
tp@0 841 firstedgeend = full(edgeextradata(ivlonglist(iv),2));
tp@0 842
tp@0 843 lastedgestart = full(edgeextradata(ivlonglist(iv),3));
tp@0 844 lastedgeend = full(edgeextradata(ivlonglist(iv),4));
tp@0 845
tp@0 846 % Calculate the cyl coordinates of all IS/S and IR/R
tp@0 847
tp@0 848 cylS = zeros(ncombs,3);
tp@0 849 edgecoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)];
tp@0 850 [cylS(:,1),cylS(:,2),cylS(:,3)] = EDB1coordtrans1(IS,edgecoords,edgenvecs(edgeshortlist(ii,1),:));
tp@0 851
tp@0 852 cylR = zeros(ncombs,3);
tp@0 853 edgecoords = [edgestartcoords(edgeshortlist(ii,Ndifforder),:);edgeendcoords(edgeshortlist(ii,Ndifforder),:)];
tp@0 854 [cylR(:,1),cylR(:,2),cylR(:,3)] = EDB1coordtrans1(IR,edgecoords,edgenvecs(edgeshortlist(ii,Ndifforder),:));
tp@0 855
tp@0 856 bc = ones(Ndifforder,2); % Check real reflfactors!!!
tp@0 857 bc = reshape(bc.',2*Ndifforder,1);
tp@0 858
tp@0 859 % Pick out the edge-to-edge coordinates for this specific
tp@0 860 % edge pair/N-let.
tp@0 861
tp@0 862 if ~isempty(reftoshortlistE),
tp@0 863 if Ndifforder >= 2
tp@0 864 index1 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,1));
tp@0 865 cylE2_r1 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 866 index2 = reftoshortlistE(edgeshortlist(ii,1),edgeshortlist(ii,2));
tp@0 867 cylE1_r2 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 868 end
tp@0 869 if Ndifforder >= 3
tp@0 870 index1 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,2));
tp@0 871 cylE3_r2 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 872 index2 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,3));
tp@0 873 cylE2_r3 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 874 end
tp@0 875 if Ndifforder >= 4
tp@0 876 index1 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,3));
tp@0 877 cylE4_r3 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 878 index2 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,4));
tp@0 879 cylE3_r4 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 880 end
tp@0 881 if Ndifforder >= 5
tp@0 882 index1 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,4));
tp@0 883 cylE5_r4 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 884 index2 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,5));
tp@0 885 cylE4_r5 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 886 end
tp@0 887 if Ndifforder >= 6
tp@0 888 index1 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,5));
tp@0 889 cylE6_r5 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 890 index2 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,6));
tp@0 891 cylE5_r6 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 892 end
tp@0 893 if Ndifforder >= 7
tp@0 894 index1 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,6));
tp@0 895 cylE7_r6 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 896 index2 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,7));
tp@0 897 cylE6_r7 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 898 end
tp@0 899 if Ndifforder >= 8
tp@0 900 index1 = reftoshortlistE(edgeshortlist(ii,8),edgeshortlist(ii,7));
tp@0 901 cylE8_r7 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)];
tp@0 902 index2 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,8));
tp@0 903 cylE7_r8 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)];
tp@0 904 end
tp@0 905
tp@0 906 else
tp@0 907 error('STRANGE TO END UP HERE??')
tp@0 908 if Ndifforder == 2
tp@0 909 cylE1_r2 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))];
tp@0 910 cylE2_r1 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))];
tp@0 911 else
tp@0 912 error(['ERROR: Geometries with a single edge can not handle difforder >= 3'])
tp@0 913 end
tp@0 914 end
tp@0 915
tp@0 916 for jj = 1:ncombs
tp@0 917
tp@0 918 if multfactors(ivlonglist(iv(jj)))
tp@0 919
tp@0 920 if Ndifforder == 2
tp@0 921
tp@0 922 cylE1_r2frac = cylE1_r2;
tp@0 923 e1length = edgelengthvec(edgeshortlist(ii,1));
tp@0 924 cylE2_r1frac = cylE2_r1;
tp@0 925 e2length = edgelengthvec(edgeshortlist(ii,2));
tp@0 926
tp@0 927 if firstedgestart(jj) ~= 0
tp@0 928 cylE1_r2frac(1,3) = cylE1_r2frac(1,3) + e1length*firstedgestart(jj);
tp@0 929 end
tp@0 930 if firstedgeend(jj) ~= 1
tp@0 931 cylE1_r2frac(2,3) = cylE1_r2frac(1,3) + e1length*(1-firstedgeend(jj));
tp@0 932 end
tp@0 933 if lastedgestart(jj) ~= 0
tp@0 934 cylE2_r1frac(1,3) = cylE2_r1frac(1,3) + e2length*lastedgestart(jj);
tp@0 935 end
tp@0 936 if lastedgeend(jj) ~= 1
tp@0 937 cylE2_r1frac(2,3) = cylE2_r1frac(1,3)+ e2length*(1-lastedgeend(jj));
tp@0 938 end
tp@0 939
tp@0 940 if midspecpattern(iv(jj))~=0
tp@0 941
tp@0 942 % For the cases with specular reflections
tp@0 943 % in-between a double diffraction we must
tp@0 944 % mirror the two edges to calculate the
tp@0 945 % relative-to-edge cylindrical coordinates.
tp@0 946
tp@0 947 nspec = sum(midspecpattern(iv(jj),:)>0);
tp@0 948
tp@0 949 % Mirror edge 2 in all the specular reflection
tp@0 950 % planes, in reversed order
tp@0 951 edgestartmirr = edgestartcoords(edgeshortlist(ii,2),:);
tp@0 952 edgeendmirr = edgeendcoords(edgeshortlist(ii,2),:);
tp@0 953 edgevector = edgeendmirr - edgestartmirr;
tp@0 954 if lastedgestart(jj) ~= 0
tp@0 955 edgestartmirr = edgestartmirr + edgevector*lastedgestart(jj);
tp@0 956 else
tp@0 957 edgestartmirr = edgestartmirr + edgevector*1e-6;
tp@0 958 end
tp@0 959 if lastedgeend(jj) ~= 1
tp@0 960 edgeendmirr = edgestartmirr + edgevector*(1-lastedgeend(jj));
tp@0 961 else
tp@0 962 edgeendmirr = edgestartmirr + edgevector*(1-1e-6);
tp@0 963 end
tp@0 964 edgerefcoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)];
tp@0 965 edgerefnvec = edgenvecs(edgeshortlist(ii,1),:);
tp@0 966
tp@0 967 % If we have a specular reflection in a plane
tp@0 968 % which is perpendicular to the edge plane, we
tp@0 969 % should nudge the mirrored edge out a bit so
tp@0 970 % that there is no 0/(2*pi) mistake
tp@0 971 if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) )
tp@0 972 vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr;
tp@0 973 edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10;
tp@0 974 vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr;
tp@0 975 edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10;
tp@0 976 end
tp@0 977 xis = [edgestartmirr;edgeendmirr];
tp@0 978 for kk = nspec:-1:1
tp@0 979 xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1);
tp@0 980 end
tp@0 981 [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec);
tp@0 982 cylE2mirr_r1 = [rstart thetastart zstart;rend thetaend zend];
tp@0 983
tp@0 984 % Mirror edge 1 in all the specular reflection
tp@0 985 % planes, in forward order
tp@0 986 edgestartmirr = edgestartcoords(edgeshortlist(ii,1),:);
tp@0 987 edgeendmirr = edgeendcoords(edgeshortlist(ii,1),:);
tp@0 988 edgevector = edgeendmirr - edgestartmirr;
tp@0 989 if firstedgestart(jj) ~= 0
tp@0 990 edgestartmirr = edgestartmirr + edgevector*firstedgestart(jj);
tp@0 991 else
tp@0 992 edgestartmirr = edgestartmirr + edgevector*1e-6;
tp@0 993 end
tp@0 994 if firstedgeend(jj) ~= 1
tp@0 995 edgeendmirr = edgestartmirr + edgevector*(1-firstedgeend(jj));
tp@0 996 else
tp@0 997 edgeendmirr = edgestartmirr + edgevector*(1-1e-6);
tp@0 998 end
tp@0 999 edgerefcoords = [edgestartcoords(edgeshortlist(ii,2),:);edgeendcoords(edgeshortlist(ii,2),:)];
tp@0 1000 edgerefnvec = edgenvecs(edgeshortlist(ii,2),:);
tp@0 1001
tp@0 1002 % Normally, when there is a specular reflection
tp@0 1003 % in-between, the diffraction path will not be
tp@0 1004 % along a plane, unless:
tp@0 1005 % 1. The same edge is involved twice, with a
tp@0 1006 % reflection in a perpendicular plane
tp@0 1007 % in-between.
tp@0 1008 % 2. See below
tp@0 1009 pathalongplanewmidspec = 0;
tp@0 1010 if edgeshortlist(ii,1) == edgeshortlist(ii,2)
tp@0 1011 if thetastart == 0 | thetastart == (2*pi-closwedangvec(edgeshortlist(ii,2)))
tp@0 1012 pathalongplanewmidspec = 1;
tp@0 1013 end
tp@0 1014 end
tp@0 1015
tp@0 1016 % If we have a specular reflection in a plane
tp@0 1017 % which is perpendicular to the edge plane, we
tp@0 1018 % should nudge the mirrored edge out a bit so
tp@0 1019 % that there is no 0/(2*pi) mistake
tp@0 1020 %
tp@0 1021 % We also have case 2 here (see above) for when
tp@0 1022 % we could have a pathalongplanewmidspec
tp@0 1023 % 2. When two different edges have a perpendicular reflection
tp@0 1024 % in-between
tp@0 1025
tp@0 1026 if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) )
tp@0 1027 vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr;
tp@0 1028 edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10;
tp@0 1029 vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr;
tp@0 1030 edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10;
tp@0 1031 pathalongplanewmidspec = 1;
tp@0 1032 end
tp@0 1033 xis = [edgestartmirr;edgeendmirr];
tp@0 1034 for kk = 1:nspec
tp@0 1035 xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1);
tp@0 1036 end
tp@0 1037 [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec);
tp@0 1038 cylE1mirr_r2 = [rstart thetastart zstart;rend thetaend zend];
tp@0 1039
tp@0 1040 [irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2mirr_r1,cylE1mirr_r2,...
tp@0 1041 nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),...
tp@0 1042 EDcalcmethod,pathalongplanewmidspec,Rstart,bc,CAIR,FSAMP);
tp@0 1043
tp@0 1044 else % .... if midspecpattern(iv(jj))~=0
tp@0 1045
tp@0 1046 [irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2_r1frac,cylE1_r2frac,...
tp@0 1047 nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),...
tp@0 1048 EDcalcmethod,pathalongplane,Rstart,bc,CAIR,FSAMP);
tp@0 1049 end % .... if midspecpattern(iv(jj))~=0
tp@0 1050
tp@0 1051 if SHOWTEXT >= 6
tp@0 1052 figure(1)
tp@0 1053 plot(irnew)
tp@0 1054 figure(2)
tp@0 1055 plot(irdiff)
tp@0 1056 pause
tp@0 1057 end
tp@0 1058 IRDIFFVEC = [IRDIFFVEC;sum(irnew)];
tp@0 1059
tp@0 1060 elseif Ndifforder == 3, % if Ndifforder == 2
tp@0 1061
tp@0 1062 for kk = 1:newndivvec(1)
tp@0 1063 if SHOWTEXT >= 5
tp@0 1064 disp([' ',int2str(kk),' of ',int2str(newndivvec(1))])
tp@0 1065 end
tp@0 1066 BIGEDGE1stvalue = (kk-0.5)./newndivvec(1);
tp@0 1067 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3];
tp@0 1068
tp@0 1069 [irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,...
tp@0 1070 nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',...
tp@0 1071 dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue);
tp@0 1072 irnewpartition = real(irnewpartition);
tp@0 1073
tp@0 1074 if kk == 1
tp@0 1075 irnew = irnewpartition;
tp@0 1076 else
tp@0 1077 lengthaddition = length(irnewpartition);
tp@0 1078 lengthaccum = length(irnew);
tp@0 1079 if lengthaddition > lengthaccum
tp@0 1080 irnew = [irnew;zeros(lengthaddition-lengthaccum,1)];
tp@0 1081 end
tp@0 1082 irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition;
tp@0 1083 end
tp@0 1084 end
tp@0 1085
tp@0 1086 if SHOWTEXT >= 6
tp@0 1087 sum(irnew)
tp@0 1088 plot(irnew)
tp@0 1089 pause
tp@0 1090 end
tp@0 1091 end
tp@0 1092 if Ndifforder == 4
tp@0 1093 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4];
tp@0 1094 elseif Ndifforder == 5
tp@0 1095 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5];
tp@0 1096 elseif Ndifforder == 6
tp@0 1097 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6];
tp@0 1098 elseif Ndifforder == 7
tp@0 1099 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7];
tp@0 1100 elseif Ndifforder == 8,
tp@0 1101 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7;cylE8_r7;cylE7_r8];
tp@0 1102 end
tp@0 1103
tp@0 1104 if Ndifforder >= 4,
tp@0 1105
tp@0 1106 for kk = 1:newndivvec(1)
tp@0 1107 if SHOWTEXT >= 5
tp@0 1108 disp([' ',int2str(kk),' of ',int2str(newndivvec(1))])
tp@0 1109 end
tp@0 1110 BIGEDGE1stvalue = (kk-0.5)./newndivvec(1);
tp@0 1111
tp@0 1112 [irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,...
tp@0 1113 nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',...
tp@0 1114 dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue);
tp@0 1115 irnewpartition = real(irnewpartition);
tp@0 1116
tp@0 1117 if kk == 1
tp@0 1118 irnew = irnewpartition;
tp@0 1119 else
tp@0 1120 lengthaddition = length(irnewpartition);
tp@0 1121 lengthaccum = length(irnew);
tp@0 1122 if lengthaddition > lengthaccum
tp@0 1123 irnew = [irnew;zeros(lengthaddition-lengthaccum,1)];
tp@0 1124 end
tp@0 1125 irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition;
tp@0 1126 end
tp@0 1127 end
tp@0 1128 end
tp@0 1129
tp@0 1130 if SHOWTEXT >= 5
tp@0 1131 varname = ['allirs',int2str(Ndifforder)];
tp@0 1132 if exist(varname) == 0
tp@0 1133 eval([varname,' = irnew;'])
tp@0 1134 else
tp@0 1135 lengthaddition = length(irnew);
tp@0 1136 eval(['lengthaccum = size(',varname,',1);'])
tp@0 1137 if lengthaddition > lengthaccum
tp@0 1138 eval([varname,' = [',varname,';zeros(lengthaddition-lengthaccum,size(',varname,',2))];'])
tp@0 1139 end
tp@0 1140 eval([varname,' = [',varname,' zeros(size(',varname,',1),1)];'])
tp@0 1141 eval([varname,'(1:length(irnew),end) = irnew;'])
tp@0 1142 end
tp@0 1143 end
tp@0 1144
tp@0 1145 % Decide whether the IR should be boosted or not
tp@0 1146
tp@0 1147 boostfactor = 1;
tp@0 1148 if souspecboost ~= 1 | recspecboost ~= 1
tp@0 1149
tp@0 1150 if prespecpresent(iv(jj)) == 1
tp@0 1151 boostfactor = souspecboost;
tp@0 1152 else
tp@0 1153 if ~isempty(Sinsideplanenumber)
tp@0 1154 if Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),1) & Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),2),
tp@0 1155 boostfactor = souspecboost;
tp@0 1156 end
tp@0 1157 end
tp@0 1158
tp@0 1159 end
tp@0 1160 if postspecpresent(iv(jj)) == 1
tp@0 1161 boostfactor = boostfactor*recspecboost;
tp@0 1162 else
tp@0 1163 if ~isempty(Rinsideplanenumber)
tp@0 1164 if Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),1) & Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),2),
tp@0 1165 boostfactor = boostfactor*recspecboost;
tp@0 1166 end
tp@0 1167 end
tp@0 1168 end
tp@0 1169
tp@0 1170 end % ... if souspecboost ~= 1 | recspecboost ~= 1
tp@0 1171
tp@0 1172 % For thin plates, we must have a boost factor!
tp@0 1173 % This is because there will be multiple equivalent
tp@0 1174 % combinations passing on the rear side of the thin plate
tp@0 1175
tp@0 1176 if all( nyvec(edgeshortlist(ii,:)) == 0.5 )
tp@0 1177 boostfactor = boostfactor*2^(Ndifforder-1);
tp@0 1178 end
tp@0 1179
tp@0 1180 if boostfactor ~= 1
tp@0 1181 irnew = irnew*boostfactor;
tp@0 1182 end
tp@0 1183
tp@0 1184 irnew = irnew*double(multfactors(ivlonglist(iv(jj))));
tp@0 1185
tp@0 1186 [ndiff,ncolsdiff] = size(irdiff);
tp@0 1187 nnew = size(irnew,1);
tp@0 1188 if saveindividualdiffirs(1) == 0
tp@0 1189 if nnew > ndiff
tp@0 1190 irdiff = [irdiff;zeros(nnew-ndiff,1)];
tp@0 1191 end
tp@0 1192 irdiff(1:nnew) = irdiff(1:nnew) + irnew;
tp@0 1193 else
tp@0 1194 if Ndifforder > ncolsdiff
tp@0 1195 irdiff = [irdiff zeros(ndiff,Ndifforder-ncolsdiff)];
tp@0 1196 ncolsdiff = ncolsdiff+1;
tp@0 1197 end
tp@0 1198 if nnew > ndiff
tp@0 1199 irdiff = [irdiff;zeros(nnew-ndiff,Ndifforder)];
tp@0 1200 end
tp@0 1201 irdiff(1:nnew,Ndifforder) = irdiff(1:nnew,Ndifforder) + irnew;
tp@0 1202 end
tp@0 1203
tp@0 1204 end % ... if multfactors(ivlonglist(iv(jj)))
tp@0 1205
tp@0 1206 end % ....for jj = 1:ncombs
tp@0 1207
tp@0 1208 else % ... if any(multfactors(ivlonglist(iv)))
tp@0 1209 if SHOWTEXT >= 4
tp@0 1210 disp([' Combination not computed because of repetition'])
tp@0 1211 end
tp@0 1212
tp@0 1213 end % ... if any(multfactors(ivlonglist(iv)))
tp@0 1214
tp@0 1215 end % .... for ii = 1: length(ivec)
tp@0 1216
tp@0 1217 end % ... if any(difforderinlist==Ndifforder)
tp@0 1218
tp@0 1219 if size(irdiff,2) < Ndifforder & saveindividualdiffirs(1) == 1
tp@0 1220 irdiff = [irdiff zeros(size(irdiff,1),Ndifforder-size(irdiff,2))];
tp@0 1221 end
tp@0 1222 end % ... for Ndifforder = 2:specorder
tp@0 1223 end % .... if userwantsdiff2 == 1
tp@0 1224
tp@0 1225 ntot = length(irtot);
tp@0 1226 ndirect = length(irdirect);
tp@0 1227 ngeom = size(irgeom,1);
tp@0 1228 ndiff = size(irdiff,1);
tp@0 1229
tp@0 1230 if ndirect > ntot
tp@0 1231 irtot = [irtot;zeros(ndirect-ntot,1)];
tp@0 1232 ntot = length(irtot);
tp@0 1233 end
tp@0 1234 irtot(1:ndirect) = irtot(1:ndirect) + irdirect;
tp@0 1235
tp@0 1236 if ngeom > ntot
tp@0 1237 irtot = [irtot;zeros(ngeom-ntot,1)];
tp@0 1238 ntot = length(irtot);
tp@0 1239 end
tp@0 1240 irtot(1:ngeom) = irtot(1:ngeom) + irgeom;
tp@0 1241
tp@0 1242 if ndiff > ntot
tp@0 1243 irtot = [irtot;zeros(ndiff-ntot,1)];
tp@0 1244 ntot = length(irtot);
tp@0 1245 end
tp@0 1246 if saveindividualdiffirs(1) == 0 | userwantsdiff2 == 0
tp@0 1247 irtot(1:ndiff) = irtot(1:ndiff) + irdiff;
tp@0 1248 else
tp@0 1249 irtot(1:ndiff) = irtot(1:ndiff) + sum(irdiff.').';
tp@0 1250 end
tp@0 1251
tp@0 1252 %-------------------------------------------------------
tp@0 1253 % Make the IRs the same length
tp@0 1254
tp@0 1255 nmax = max([ length(irtot) size(irgeom,1) size(irdiff,1) length(irdirect)]);
tp@0 1256
tp@0 1257 if length(irtot) < nmax
tp@0 1258 irtot = [irtot;zeros(nmax-length(irtot),1)];
tp@0 1259 end
tp@0 1260 if length(irdirect) < nmax
tp@0 1261 irdirect = [irdirect;zeros(nmax-length(irdirect),1)];
tp@0 1262 end
tp@0 1263 if length(irgeom) < nmax
tp@0 1264 irgeom = [irgeom;zeros(nmax-size(irgeom,1),size(irgeom,2))];
tp@0 1265 end
tp@0 1266 if length(irdiff) < nmax
tp@0 1267 irdiff = [irdiff;zeros(nmax-size(irdiff,1),size(irdiff,2))];
tp@0 1268 end
tp@0 1269
tp@0 1270 %-------------------------------------------------------
tp@0 1271 % Save the IRs
tp@0 1272
tp@0 1273 irtot = sparse(irtot);
tp@0 1274 irgeom = sparse(irgeom);
tp@0 1275 irdirect = sparse(irdirect);
tp@0 1276 irdiff = sparse(irdiff);
tp@0 1277 Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod'];
tp@0 1278
tp@0 1279 if length(saveindividualdiffirs) > 1
tp@0 1280 if saveindividualdiffirs(2) == 1
tp@0 1281 Varlist = [Varlist,' alldiffirs alldiffirsextradata'];
tp@0 1282 end
tp@0 1283 end
tp@0 1284
tp@0 1285 if SHOWTEXT >= 5
tp@0 1286 if specorder >= 2
tp@0 1287 if exist('allirs2') == 0
tp@0 1288 allirs2 = [];
tp@0 1289 end
tp@0 1290 Varlist = [Varlist,' allirs2'];
tp@0 1291 end
tp@0 1292 if specorder >= 3
tp@0 1293 if exist('allirs3') == 0
tp@0 1294 allirs3 = [];
tp@0 1295 end
tp@0 1296 Varlist = [Varlist,' allirs3'];
tp@0 1297 end
tp@0 1298 if specorder >= 4
tp@0 1299 if exist('allirs4') == 0
tp@0 1300 allirs4 = [];
tp@0 1301 end
tp@0 1302 Varlist = [Varlist,' allirs4'];
tp@0 1303 end
tp@0 1304 if specorder >= 5
tp@0 1305 if exist('allirs5') == 0
tp@0 1306 allirs5 = [];
tp@0 1307 end
tp@0 1308 Varlist = [Varlist,' allirs5'];
tp@0 1309 end
tp@0 1310 if specorder >= 6
tp@0 1311 if exist('allirs6') == 0
tp@0 1312 allirs6 = [];
tp@0 1313 end
tp@0 1314 Varlist = [Varlist,' allirs6'];
tp@0 1315 end
tp@0 1316 if specorder >= 7
tp@0 1317 if exist('allirs7') == 0
tp@0 1318 allirs7 = [];
tp@0 1319 end
tp@0 1320 Varlist = [Varlist,' allirs7'];
tp@0 1321 end
tp@0 1322 if specorder >= 8
tp@0 1323 if exist('allirs8') == 0
tp@0 1324 allirs8 = [];
tp@0 1325 end
tp@0 1326 Varlist = [Varlist,' allirs8'];
tp@0 1327 end
tp@0 1328 end
tp@0 1329
tp@0 1330 eval(['save ',edirfile,Varlist])