Mercurial > hg > human-echolocation

annotate private/EDB1findISEStree.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 ISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfromS,vispartedgesfromS,nedgesubs)
tp@0 2 % EDB1findISEStree - Constructs a list ("an ISES tree") of all possible specular-diffraction combinations.
tp@0 3 % EDB1findISEStree constructs a list ("an ISES tree") of all possible specular-diffraction combinations that
tp@0 4 % are possible for the given source. This list is based on the visibility
tp@0 5 % matrices from source to planes/edges and between planes/edges. Visibility
tp@0 6 % checks are employed wherever possibble. NB! Visibility checks that depend
tp@0 7 % on the receiver position are not used at this stage.
tp@0 8 % Planes are numbered 1-nplanes, and edges have numbers
tp@0 9 % [nplanes+1:nplanes+nedges].
tp@0 10 %
tp@0 11 % Input parameters:
tp@0 12 % eddatafile See a description in EDB1edgeo, EDB1srgeo and EDB1mainISES
tp@0 13 % S
tp@0 14 % isou
tp@0 15 % specorder
tp@0 16 % difforder
tp@0 17 % visplanesfromS
tp@0 18 % vispartedgesfromS
tp@0 19 % nedgesubs
tp@0 20 %
tp@0 21 % GLOBAL parameters
tp@0 22 % SHOWTEXT JJ JJnumbofchars See EDB1mainISES
tp@0 23 % SUPPRESSFILES If this global parameter has the value 1
tp@0 24 % then the results from this function will be returned in a struct
tp@0 25 % rather than in a file.
tp@0 26 %
tp@0 27 % Output parameters:
tp@0 28 % ISEStreefile The name of the output file which contains
tp@0 29 % the output data below.
tp@0 30 %
tp@0 31 % Global output data:
tp@0 32 % POTENTIALISES A matrix, [npossiblecombs,specorder], of all the
tp@0 33 % possible specular-diffraction combos that are possible for the
tp@0 34 % given source. The combos are denoted by
tp@0 35 % plane numbers and edge numbers, but edge
tp@0 36 % numbers have a constant number (=nplanes)
tp@0 37 % added to them, i.e., if there are 20 planes
tp@0 38 % in the model, edge number 1 will be denoted
tp@0 39 % by number 21.
tp@0 40 % ORIGINSFROM A list, [npossiblecombs,1], where the value
tp@0 41 % in row N states which other row in ISCOORDS that
tp@0 42 % the image source in row N originates from.
tp@0 43 % ISCOORDS A matrix, [npossiblecombs,3], containing
tp@0 44 % the image source coordinates, where this is applicable
tp@0 45 % i.e., where the combo starts with a specular reflection.
tp@0 46 % ISESVISIBILITY A list, [npossiblecombs,1], containing the visibility of
tp@0 47 % the first edge in a sequence, seen from the source.
tp@0 48 % IVNSPECMATRIX A matrix, [nspeccombs,specorder], where each column contains
tp@0 49 % the row numbers in POTENTIALISES that contain combos with
tp@0 50 % 1 or 2 or... or "specorder" purely specular reflections.
tp@0 51 % REFLORDER A list, [npossiblecombs,1], containing the
tp@0 52 % order of reflection for each row of POTENTIALISES
tp@0 53 % IVNDIFFMATRIX A matrix, [ndiffcombs,specorder], where each column contains
tp@0 54 % the row numbers in POTENTIALISES that contain combos with
tp@0 55 % 1 or 2 or... or "specorder" diffractions.
tp@0 56 % Output data in the ISEStreefile:
tp@0 57 % lengthNspecmatrix A list, [1,specorder], with the number of
tp@0 58 % entries in each column of IVNSPECMATRIX.
tp@0 59 % lengthNdiffmatrix A list, [1,specorder], with the number of
tp@0 60 % entries in each column of IVNDIFFMATRIX.
tp@0 61 % singlediffcol A list, [nfirstorderdiff,1], containing the column number
tp@0 62 % that the first-order diffracted component can be found in (in
tp@0 63 % POTENTIALISES).
tp@0 64 % startindicessinglediff A list, [specorder,1], where the value in position N contains
tp@0 65 % which row in IVNDIFFMATRIX(:,1) that has the first combination
tp@0 66 % of order N and with one diffraction component.
tp@0 67 % endindicessinglediff A list, [specorder,1], where the value in position N contains
tp@0 68 % which row in IVNDIFFMATRIX(:,1) that has the last combination
tp@0 69 % of order N and with one diffraction
tp@0 70 % ndecimaldivider See below.
tp@0 71 % PointertoIRcombs A sparse list which for position N contains a row number
tp@0 72 % (in POTENTIALISES) that has a combination
tp@0 73 % that ends with a specular reflection in
tp@0 74 % plane N, after a diffraction. A row with one of the specular
tp@0 75 % reflections M,N after a diffraction can be
tp@0 76 % found in position M*ndecimaldivider+N etc
tp@0 77 % for higher orders.
tp@0 78 % IRoriginsfrom A list, [npossiblecombs,1], where the value in position N
tp@0 79 % states which other row number (in POTENTIALISES) that the
tp@0 80 % image receiver in row N origins from.
tp@0 81 %
tp@0 82 % Uses functions EDB1strpend EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths
tp@0 83 %
tp@0 84 % ----------------------------------------------------------------------------------------------
tp@0 85 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 86 %
tp@0 87 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 88 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 89 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 90 %
tp@0 91 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 92 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 93 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 94 %
tp@0 95 % You should have received a copy of the GNU General Public License along with the
tp@0 96 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 97 % ----------------------------------------------------------------------------------------------
tp@0 98 % Peter Svensson (svensson@iet.ntnu.no) 20100210
tp@0 99 %
tp@0 100 % ISEStreefile = EDB1findISEStree(eddatafile,S,isou,specorder,difforder,visplanesfroms,vispartedgesfroms,nedgesubs);
tp@0 101
tp@0 102 global SHOWTEXT JJ JJnumbofchars
tp@0 103 global POTENTIALISES ISCOORDS IVNDIFFMATRIX
tp@0 104 global IVNSPECMATRIX ORIGINSFROM ISESVISIBILITY REFLORDER SUPPRESSFILES
tp@0 105
tp@0 106 GEOMACC = 1e-10;
tp@0 107
tp@0 108 eval(['load ',eddatafile])
tp@0 109 clear cornerinfrontofplane
tp@0 110
tp@0 111 [filepath,filestem,fileext] = fileparts(eddatafile);
tp@0 112 ISEStreefile = [[filepath,filesep],EDB1strpend(filestem,'_eddata'),'_',int2str(isou),'_ISEStree.mat'];
tp@0 113
tp@0 114 ncorners = size(corners,1);
tp@0 115 nplanes = length(planeisthin);
tp@0 116 nedges = length(closwedangvec);
tp@0 117 nhighval = nplanes + nedges;
tp@0 118
tp@0 119 onesvec1 = ones(1,nedgesubs);
tp@0 120 obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0 121 onesvec3 = ones(1,3);
tp@0 122
tp@0 123 doconetrace = 0;
tp@0 124 if SHOWTEXT >= 3
tp@0 125 disp('WARNING: doconetrace is switched off')
tp@0 126 end
tp@0 127
tp@0 128 %--------------------------------------------------------------------------
tp@0 129 % We modify edgeseesplane because in EDB1edgeo, edgeseesplane was set to 1
tp@0 130 % for totabs planes (this was necessary for EDB1srgeo).
tp@0 131
tp@0 132 % Set all edges belonging to planes with reflfactors = 0 to edgeseesplane = 0.
tp@0 133
tp@0 134 listofabsorbingplanes = find(reflfactors == 0);
tp@0 135
tp@0 136 if ~isempty(listofabsorbingplanes)
tp@0 137 edgeseesplane(listofabsorbingplanes,:) = zeros(length(listofabsorbingplanes),nedges);
tp@0 138 end
tp@0 139
tp@0 140 % If we should do cone tracing, then we need plane midpoints and plane
tp@0 141 % sizes.
tp@0 142
tp@0 143 if doconetrace == 1
tp@0 144 planemidpoints = zeros(nplanes,3);
tp@0 145 maxdisttocorner = zeros(nplanes,1);
tp@0 146
tp@0 147 iv4planes = find(ncornersperplanevec==4);
tp@0 148 if any(iv4planes)
tp@0 149 xvalues = corners(planecorners(iv4planes,1:4).',1);
tp@0 150 xvalues = reshape(xvalues,4,length(xvalues)/4);
tp@0 151 yvalues = corners(planecorners(iv4planes,1:4).',2);
tp@0 152 yvalues = reshape(yvalues,4,length(yvalues)/4);
tp@0 153 zvalues = corners(planecorners(iv4planes,1:4).',3);
tp@0 154 zvalues = reshape(zvalues,4,length(zvalues)/4);
tp@0 155 planemidpoints(iv4planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0 156
tp@0 157 dist1 = sqrt(sum((corners(planecorners(iv4planes,1),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0 158 dist2 = sqrt(sum((corners(planecorners(iv4planes,2),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0 159 dist3 = sqrt(sum((corners(planecorners(iv4planes,3),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0 160 dist4 = sqrt(sum((corners(planecorners(iv4planes,4),:) - planemidpoints(iv4planes,:)).'.^2)).';
tp@0 161 maxdisttocorner(iv4planes) = max([dist1 dist2 dist3 dist4].');
tp@0 162
tp@0 163 end
tp@0 164
tp@0 165 iv3planes = find(ncornersperplanevec==3);
tp@0 166 if any(iv3planes)
tp@0 167 xvalues = corners(planecorners(iv3planes,1:3).',1);
tp@0 168 xvalues = reshape(xvalues,3,length(xvalues)/3);
tp@0 169 yvalues = corners(planecorners(iv3planes,1:3).',2);
tp@0 170 yvalues = reshape(yvalues,3,length(yvalues)/3);
tp@0 171 zvalues = corners(planecorners(iv3planes,1:3).',3);
tp@0 172 zvalues = reshape(zvalues,3,length(zvalues)/3);
tp@0 173 planemidpoints(iv3planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0 174
tp@0 175 dist1 = sqrt(sum((corners(planecorners(iv3planes,1),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0 176 dist2 = sqrt(sum((corners(planecorners(iv3planes,2),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0 177 dist3 = sqrt(sum((corners(planecorners(iv3planes,3),:) - planemidpoints(iv3planes,:)).'.^2)).';
tp@0 178 maxdisttocorner(iv3planes) = max([dist1 dist2 dist3].');
tp@0 179
tp@0 180 end
tp@0 181
tp@0 182 iv5planes = find(ncornersperplanevec==5);
tp@0 183 if any(iv5planes)
tp@0 184 xvalues = corners(planecorners(iv5planes,1:5).',1);
tp@0 185 xvalues = reshape(xvalues,5,length(xvalues)/5);
tp@0 186 yvalues = corners(planecorners(iv5planes,1:5).',2);
tp@0 187 yvalues = reshape(yvalues,5,length(yvalues)/5);
tp@0 188 zvalues = corners(planecorners(iv5planes,1:5).',3);
tp@0 189 zvalues = reshape(zvalues,5,length(zvalues)/5);
tp@0 190 planemidpoints(iv5planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0 191
tp@0 192 dist1 = sqrt(sum((corners(planecorners(iv5planes,1),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0 193 dist2 = sqrt(sum((corners(planecorners(iv5planes,2),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0 194 dist3 = sqrt(sum((corners(planecorners(iv5planes,3),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0 195 dist4 = sqrt(sum((corners(planecorners(iv5planes,4),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0 196 dist5 = sqrt(sum((corners(planecorners(iv5planes,5),:) - planemidpoints(iv5planes,:)).'.^2)).';
tp@0 197 maxdisttocorner(iv5planes) = max([dist1 dist2 dist3 dist4 dist5].');
tp@0 198
tp@0 199 end
tp@0 200
tp@0 201 iv6planes = find(ncornersperplanevec==6);
tp@0 202 if any(iv6planes)
tp@0 203 xvalues = corners(planecorners(iv6planes,1:6).',1);
tp@0 204 xvalues = reshape(xvalues,6,length(xvalues)/6);
tp@0 205 yvalues = corners(planecorners(iv6planes,1:6).',2);
tp@0 206 yvalues = reshape(yvalues,6,length(yvalues)/6);
tp@0 207 zvalues = corners(planecorners(iv6planes,1:6).',3);
tp@0 208 zvalues = reshape(zvalues,6,length(zvalues)/6);
tp@0 209 planemidpoints(iv6planes,:) = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
tp@0 210
tp@0 211 dist1 = sqrt(sum((corners(planecorners(iv6planes,1),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 212 dist2 = sqrt(sum((corners(planecorners(iv6planes,2),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 213 dist3 = sqrt(sum((corners(planecorners(iv6planes,3),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 214 dist4 = sqrt(sum((corners(planecorners(iv6planes,4),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 215 dist5 = sqrt(sum((corners(planecorners(iv6planes,5),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 216 dist6 = sqrt(sum((corners(planecorners(iv6planes,6),:) - planemidpoints(iv6planes,:)).'.^2)).';
tp@0 217 maxdisttocorner(iv6planes) = max([dist1 dist2 dist3 dist4 dist5 dist6].');
tp@0 218
tp@0 219 end
tp@0 220
tp@0 221 end
tp@0 222
tp@0 223 %--------------------------------------------------------------------------
tp@0 224 % Combine planeseesplane with edgeseesplane
tp@0 225 % and visplanesfroms with vispartedgesfroms
tp@0 226 % so that edges appear as planes
tp@0 227
tp@0 228 if difforder > 0
tp@0 229 visPLANESfroms = [visplanesfromS;2*sign(double(vispartedgesfromS))];
tp@0 230 % source sees planes that have visplanesfroms:
tp@0 231 % 2 => in front of reflective planes
tp@0 232 % 4 => inside plane which is reflective
tp@0 233 % old version: sousees = (visPLANESfroms==1 | visPLANESfroms==-2);
tp@0 234 sousees = (visPLANESfroms==2 | visPLANESfroms==4);
tp@0 235
tp@0 236 if exist('edgeseespartialedge') ~= 1
tp@0 237 edgeseesedge = int8(zeros(nedges,nedges));
tp@0 238 else
tp@0 239 if ~isempty(edgeseespartialedge)
tp@0 240
tp@0 241 edgeseesedge = int8(full(abs(double(edgeseespartialedge))>0));
tp@0 242
tp@0 243 else
tp@0 244 edgeseesedge = int8(zeros(nedges,nedges));
tp@0 245 end
tp@0 246 end
tp@0 247
tp@0 248 PLANEseesPLANE = [[planeseesplane edgeseesplane];[edgeseesplane.' edgeseesedge]];
tp@0 249
tp@0 250 clear edgeseesplane edgeseesedge
tp@0 251 else
tp@0 252 visPLANESfroms = visplanesfromS;
tp@0 253 PLANEseesPLANE = planeseesplane;
tp@0 254 sousees = (visPLANESfroms==2 | visPLANESfroms==4);
tp@0 255
tp@0 256 clear visplanesfromS
tp@0 257 end
tp@0 258
tp@0 259 startindices = zeros(specorder,1);
tp@0 260 endindices = zeros(specorder,1);
tp@0 261
tp@0 262 %--------------------------------------------------------------------------
tp@0 263 % Construct a list of which planes a plane sees so that a search isn't
tp@0 264 % needed later.
tp@0 265
tp@0 266 nPLANES = size(PLANEseesPLANE,1);
tp@0 267 listofvisPLANES = zeros(nPLANES,nPLANES-1);
tp@0 268 nvisPLANES = zeros(nPLANES,1);
tp@0 269 for ii = 1:nPLANES
tp@0 270 visPLANES = find(PLANEseesPLANE(ii,:)==1);
tp@0 271 nvisPLANES(ii) = length(visPLANES);
tp@0 272 listofvisPLANES(ii,1:nvisPLANES(ii)) = visPLANES;
tp@0 273 end
tp@0 274
tp@0 275 %##################################################################
tp@0 276 %##################################################################
tp@0 277 %##################################################################
tp@0 278 %
tp@0 279 % First order
tp@0 280 %
tp@0 281 %------------------------------------------------------------------
tp@0 282
tp@0 283 startindices(1) = 1;
tp@0 284
tp@0 285 possiblefirsts = find( sousees );
tp@0 286 npossiblefirsts = length(possiblefirsts);
tp@0 287
tp@0 288 if nhighval < 255
tp@0 289 POTENTIALISES = uint8( [[possiblefirsts ] zeros(npossiblefirsts,specorder-1)] );
tp@0 290 else
tp@0 291 POTENTIALISES = uint16( [[possiblefirsts ] zeros(npossiblefirsts,specorder-1)] );
tp@0 292 end
tp@0 293
tp@0 294 ORIGINSFROM = uint32(zeros(npossiblefirsts,1));
tp@0 295 if nedgesubs < 8
tp@0 296 ISESVISIBILITY = uint8(zeros(npossiblefirsts,1));
tp@0 297 elseif nedgesubs < 16
tp@0 298 ISESVISIBILITY = uint16(zeros(npossiblefirsts,1));
tp@0 299 else
tp@0 300 ISESVISIBILITY = uint32(zeros(npossiblefirsts,1));
tp@0 301 end
tp@0 302 iv = find(possiblefirsts>nplanes);
tp@0 303 if ~isempty(iv)
tp@0 304 ISESVISIBILITY(iv) = vispartedgesfromS(possiblefirsts(iv)-nplanes);
tp@0 305 end
tp@0 306
tp@0 307 endindices(1) = npossiblefirsts;
tp@0 308
tp@0 309 % Compute the IS coords
tp@0 310
tp@0 311 ivdiff = [startindices(1):endindices(1)];
tp@0 312 ivspec = find(POTENTIALISES(ivdiff,1)<=nplanes);
tp@0 313 ivdiff(ivspec) = [];
tp@0 314
tp@0 315 ISCOORDS = zeros(npossiblefirsts,3);
tp@0 316 ISCOORDS(ivspec,:) = EDB1findis(S,POTENTIALISES(ivspec,1),planeeqs,1,onesvec3);
tp@0 317
tp@0 318 %##################################################################
tp@0 319 %##################################################################
tp@0 320 %##################################################################
tp@0 321 %
tp@0 322 % Second order
tp@0 323 %
tp@0 324 %---------------------------------------------------------------------------------------------
tp@0 325
tp@0 326 startindices(2) = startindices(1) + length(possiblefirsts);
tp@0 327
tp@0 328 if specorder > 1
tp@0 329 if SHOWTEXT >= 3
tp@0 330 disp([' Order number 2'])
tp@0 331 end
tp@0 332
tp@0 333 if nhighval < 255
tp@0 334 startplanes = uint8(possiblefirsts);
tp@0 335 else
tp@0 336 startplanes = uint16(possiblefirsts);
tp@0 337 end
tp@0 338
tp@0 339 addtocomb = startplanes;
tp@0 340 nadds = size(addtocomb,1);
tp@0 341 if nadds > 0
tp@0 342
tp@0 343 maxnumberofvispla = max(sum(PLANEseesPLANE(:,startplanes)==1));
tp@0 344
tp@0 345 addtocomb = (reshape(addtocomb(:,ones(1,maxnumberofvispla)).',length(addtocomb)*maxnumberofvispla,1));
tp@0 346 addtocomb = [addtocomb zeros(size(addtocomb))];
tp@0 347
tp@0 348 addtoISESVISIBILITY = reshape(ISESVISIBILITY(:,ones(1,maxnumberofvispla)).',nadds*maxnumberofvispla,1);
tp@0 349
tp@0 350 startpos = [[1:maxnumberofvispla:length(addtocomb)] length(addtocomb)+1 ].';
tp@0 351
tp@0 352 for ii = 1:length(startpos)-1
tp@0 353 if nvisPLANES(addtocomb(startpos(ii))) > 0
tp@0 354 possibleplanes = listofvisPLANES(addtocomb(startpos(ii)),1:nvisPLANES(addtocomb(startpos(ii)))).';
tp@0 355 nposs = length(possibleplanes);
tp@0 356 addtocomb( startpos(ii):startpos(ii)+nposs-1,2) = possibleplanes;
tp@0 357 end
tp@0 358
tp@0 359 end
tp@0 360
tp@0 361 addtooriginsfrom = uint32([1:startindices(2)-1].');
tp@0 362 addtooriginsfrom = addtooriginsfrom(:,ones(1,maxnumberofvispla));
tp@0 363 addtooriginsfrom = reshape(addtooriginsfrom.',maxnumberofvispla*(startindices(2)-1),1);
tp@0 364
tp@0 365 cutout = find(addtocomb(:,2)==0);
tp@0 366 addtocomb(cutout,:) = [];
tp@0 367 addtooriginsfrom(cutout) = [];
tp@0 368 addtoISESVISIBILITY(cutout) = [];
tp@0 369 clear cutout
tp@0 370 nadds = size(addtocomb,1);
tp@0 371
tp@0 372 if nadds > 0
tp@0 373
tp@0 374 %--------------------------------------------------------------------------
tp@0 375 % Try to get rid of some combinations that we know can not be propagated
tp@0 376
tp@0 377 % Find those combinations of specular-specular where the IS is behind the second reflecting plane
tp@0 378 % because they can then never ever be propagated, or valid
tp@0 379
tp@0 380 if difforder > 0
tp@0 381 ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0 382 else
tp@0 383 ivss = uint32([1:nadds].');
tp@0 384 end
tp@0 385
tp@0 386 imsoucheck = dot((ISCOORDS(addtooriginsfrom(ivss),1:3) - corners(planecorners(addtocomb(ivss,2),1),:)).',planenvecs(addtocomb(ivss,2),:).').';
tp@0 387 imsounotOK = find(imsoucheck < GEOMACC);
tp@0 388 addtocomb(ivss(imsounotOK),:) = [];
tp@0 389 addtooriginsfrom(ivss(imsounotOK),:) = [];
tp@0 390 addtoISESVISIBILITY(ivss(imsounotOK)) = [];
tp@0 391
tp@0 392 nadds = size(addtocomb,1);
tp@0 393
tp@0 394 if nadds > 0 & doconetrace == 1
tp@0 395
tp@0 396 if difforder > 0
tp@0 397 ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0 398 else
tp@0 399 ivss = uint32([1:nadds].');
tp@0 400 end
tp@0 401
tp@0 402 beamdirection1 = planemidpoints(addtocomb(ivss,1),:) - ISCOORDS(addtocomb(ivss,1),:);
tp@0 403 beamlength = sqrt( sum( beamdirection1.'.^2 ) ).';
tp@0 404 beamdirection1 = beamdirection1./beamlength(:,ones(1,3));
tp@0 405 maxprojradius = maxdisttocorner(addtocomb(ivss,1));
tp@0 406 % The maxdisttocorner below should be multiplied
tp@0 407 % by sine of the angle between the plane normal vector and the
tp@0 408 % beam direction.
tp@0 409 iv = find(beamlength<maxprojradius);
tp@0 410 if ~isempty(iv)
tp@0 411 beamlength(iv) = beamlength(iv)*0+ eps*1e6;
tp@0 412 end
tp@0 413 beamangle1 = atan(maxprojradius./beamlength);
tp@0 414
tp@0 415 % Now we calculate the beam-widths to the second reflection
tp@0 416 % plane, seen from the first-order image source. If this
tp@0 417 % beamwidth is outside the previously calculated beam of
tp@0 418 % reflection plane 1, then plane 2 can never be seen
tp@0 419 % through plane 1.
tp@0 420
tp@0 421 beamdirection2 = planemidpoints(addtocomb(ivss,2),:) - ISCOORDS(addtocomb(ivss,1),:);
tp@0 422 beamlength = sqrt( sum( beamdirection2.'.^2 ) ).';
tp@0 423 beamdirection2 = beamdirection2./beamlength(:,ones(1,3));
tp@0 424 maxprojradius = maxdisttocorner(addtocomb(ivss,2));
tp@0 425 % The maxdisttocorner below should be multiplied
tp@0 426 % by sine of the angle between the plane normal vector and the
tp@0 427 % beam direction.
tp@0 428 iv = find(beamlength<maxprojradius);
tp@0 429 if ~isempty(iv)
tp@0 430 beamlength(iv) = beamlength(iv)*0+ eps*10;
tp@0 431 end
tp@0 432 beamangle2 = atan(maxprojradius./beamlength);
tp@0 433 anglebetweenbeams = acos(sum( (beamdirection1.*beamdirection2).' ).');
tp@0 434 ivcutout = find( (beamangle1 + beamangle2 < anglebetweenbeams) & (beamangle1 <1.56) & (beamangle2 < 1.56));
tp@0 435 if ~isempty(ivcutout)
tp@0 436 addtocomb(ivss(ivcutout),:) = [];
tp@0 437 addtooriginsfrom(ivss(ivcutout),:) = [];
tp@0 438 addtoISESVISIBILITY(ivss(ivcutout)) = [];
tp@0 439 nadds = size(addtocomb,1);
tp@0 440 end
tp@0 441
tp@0 442 end
tp@0 443
tp@0 444 % Combinations of diffraction-specular don't need to be checked because
tp@0 445 % 'edgeseesplane' should have taken care of this.
tp@0 446 %%%%ivds = find(addtocomb(:,1)>nplanes & addtocomb(:,2)<=nplanes);
tp@0 447
tp@0 448 %--------------------------------------------------------------------------
tp@0 449 % Combinations of specular-diffraction: check if the IS is behind both of
tp@0 450 % the two planes that the reflecting edge connects to.
tp@0 451 %
tp@0 452 % Bug found 20080711 PS: The visibility test must be split into two, depending on the wedge angle!!
tp@0 453 % If the open wedge angle > pi, then it is correct to check if the IS behind both of the planes.
tp@0 454 % If the open wedge angle < pi, then it should be checked if the IS behind one of the planes!!
tp@0 455
tp@0 456 if difforder > 0 & nadds > 0
tp@0 457
tp@0 458 % First we check the sd combinations where the open wedge angle > pi
tp@0 459 % For those cases, we can exclude combinations where the IS
tp@0 460 % is behind both planes
tp@0 461
tp@0 462 ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0 463 ivsd1 = ivsd( find( closwedangvec( double(addtocomb(ivsd,2))-nplanes ) < pi ));
tp@0 464
tp@0 465 if ~isempty(ivsd1)
tp@0 466 imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,2))-nplanes,1) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,2))-nplanes,1) ,:).').';
tp@0 467 imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,2))-nplanes,2) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,2))-nplanes,2) ,:).').';
tp@0 468 GEOMACC = 1e-10;
tp@0 469 imsounotOK = find(imsoucheck1 < GEOMACC & imsoucheck2 < GEOMACC);
tp@0 470 clear imsoucheck1 imsoucheck2
tp@0 471
tp@0 472 addtocomb(ivsd1(imsounotOK),:) = [];
tp@0 473 addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0 474 addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0 475 end
tp@0 476
tp@0 477 % Second we check the sd combinations where the open wedge angle < pi
tp@0 478 % For those cases, we can exclude combinations where the IS
tp@0 479 % is behind one of the two planes
tp@0 480
tp@0 481 ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0 482 if ~isempty(ivsd)
tp@0 483 ivsd1 = ivsd( find( closwedangvec( double(addtocomb(ivsd,2))-nplanes ) > pi ));
tp@0 484
tp@0 485 if ~isempty(ivsd1)
tp@0 486 imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,2))-nplanes,1) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,2))-nplanes,1) ,:).').';
tp@0 487 imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,2))-nplanes,2) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,2))-nplanes,2) ,:).').';
tp@0 488 GEOMACC = 1e-10;
tp@0 489 imsounotOK = find(imsoucheck1 < GEOMACC | imsoucheck2 < GEOMACC);
tp@0 490 clear imsoucheck1 imsoucheck2
tp@0 491
tp@0 492 addtocomb(ivsd1(imsounotOK),:) = [];
tp@0 493 addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0 494 addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0 495 end
tp@0 496 end
tp@0 497
tp@0 498 nadds = size(addtocomb,1);
tp@0 499
tp@0 500 ivsd = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes));
tp@0 501
tp@0 502 if ~isempty(ivsd)
tp@0 503
tp@0 504 %----------------------------------------------------------------------
tp@0 505 % Combinations of specular-diffraction that are not visible at
tp@0 506 % all should be removed and then they are not propagated either.
tp@0 507
tp@0 508 if nedges < 255
tp@0 509 possibleedges = uint8(double(addtocomb(ivsd,2)) - nplanes);
tp@0 510 else
tp@0 511 possibleedges = uint16(double(addtocomb(ivsd,2)) - nplanes);
tp@0 512 end
tp@0 513
tp@0 514 possiblecombs = addtocomb(ivsd,1);
tp@0 515 reftoISCOORDS = addtooriginsfrom(ivsd);
tp@0 516
tp@0 517 % Expand to take the edge subdivisions into account
tp@0 518
tp@0 519 nposs = length(ivsd);
tp@0 520 nposs = nposs*nedgesubs; % We treat all the little edge subdivisions as separate edges
tp@0 521
tp@0 522 expandedposscombs = possiblecombs(:,onesvec1);
tp@0 523 clear possiblecombs
tp@0 524 expandedposscombs = reshape(expandedposscombs.',nposs,1);
tp@0 525
tp@0 526 expandedreftoISCOORDS = reftoISCOORDS(:,onesvec1);
tp@0 527 expandedreftoISCOORDS = reshape(expandedreftoISCOORDS.',nposs,1);
tp@0 528 expandedpossedges = possibleedges(:,onesvec1);
tp@0 529 expandedpossedges = reshape(expandedpossedges.',nposs,1);
tp@0 530 expandedivsd = ivsd(:,onesvec1);
tp@0 531 expandedivsd = reshape(expandedivsd.',nposs,1);
tp@0 532
tp@0 533 if SHOWTEXT >= 3
tp@0 534 disp([' ',int2str(nposs),' IS+edge segments found initially '])
tp@0 535 end
tp@0 536
tp@0 537 if nposs > 0
tp@0 538
tp@0 539 fromcoords = ISCOORDS(expandedreftoISCOORDS,:);
tp@0 540 [tocoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0 541 clear possibleedges
tp@0 542
tp@0 543 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,1),4),planenvecs(expandedposscombs(:,1),:),minvals(expandedposscombs(:,1),:),...
tp@0 544 maxvals(expandedposscombs(:,1),:),planecorners(expandedposscombs(:,1),:),corners,ncornersperplanevec(expandedposscombs(:,1)));
tp@0 545 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 546 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 547 disp(' handled correctly yet.')
tp@0 548 end
tp@0 549
tp@0 550 eval(['reflpoints',JJ(1,1),' = reflpoints;'])
tp@0 551
tp@0 552 expandedivsd = expandedivsd(hitplanes);
tp@0 553 expandedposscombs = expandedposscombs(hitplanes,:);
tp@0 554 expandedreftoISCOORDS = expandedreftoISCOORDS(hitplanes);
tp@0 555 expandedpossedges = expandedpossedges(hitplanes);
tp@0 556 edgeweightlist = edgeweightlist(hitplanes);
tp@0 557 toedgecoords = tocoords(hitplanes,:);
tp@0 558
tp@0 559 nposs = length(expandedivsd);
tp@0 560
tp@0 561 end
tp@0 562
tp@0 563 if SHOWTEXT >= 3
tp@0 564 disp([' ',int2str(nposs),' IS+edge segments visible '])
tp@0 565 end
tp@0 566
tp@0 567 if obstructtestneeded & nposs > 0
tp@0 568
tp@0 569 for jj = 1:2
tp@0 570 if nposs > 0
tp@0 571 if jj==1
tp@0 572 fromcoords = S;
tp@0 573 startplanes = [];
tp@0 574 else
tp@0 575 startplanes = expandedposscombs(:,jj-1);
tp@0 576 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0 577 end
tp@0 578 if jj == 2
tp@0 579 tocoords = toedgecoords;
tp@0 580 endplanes = [planesatedge(expandedpossedges,1) planesatedge(expandedpossedges,2)];
tp@0 581 else
tp@0 582 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
tp@0 583 endplanes = expandedposscombs(:,jj);
tp@0 584 end
tp@0 585
tp@0 586 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 587 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 588
tp@0 589 if nobstructions > 0
tp@0 590 expandedivsd = expandedivsd(nonobstructedpaths);
tp@0 591 expandedposscombs = expandedposscombs(nonobstructedpaths,:);
tp@0 592 expandedreftoISCOORDS = expandedreftoISCOORDS(nonobstructedpaths);
tp@0 593 expandedpossedges = expandedpossedges(nonobstructedpaths);
tp@0 594 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 595 toedgecoords = tocoords(nonobstructedpaths,:);
tp@0 596 nposs = length(expandedivsd);
tp@0 597
tp@0 598 for kk = 1:1, %norder
tp@0 599 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])
tp@0 600 end
tp@0 601
tp@0 602 end
tp@0 603 end
tp@0 604
tp@0 605 end
tp@0 606 end
tp@0 607
tp@0 608 if SHOWTEXT >= 3
tp@0 609 disp([' ',int2str(nposs),' IS+edge segments survived the obstruction test'])
tp@0 610 end
tp@0 611
tp@0 612 % There are repetitions of each combination since each edge segment has
tp@0 613 % been treated separately. Pack them together now
tp@0 614
tp@0 615 test = [expandedposscombs expandedpossedges];
tp@0 616 ncombs = length(expandedpossedges);
tp@0 617 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 618 ivremove = find(dtest==1);
tp@0 619
tp@0 620 while ~isempty(ivremove) & ~isempty(edgeweightlist),
tp@0 621 edgeweightlist(ivremove+1) = double(edgeweightlist(ivremove+1)) + double(edgeweightlist(ivremove));
tp@0 622 edgeweightlist(ivremove) = [];
tp@0 623 expandedpossedges(ivremove) = [];
tp@0 624 expandedposscombs(ivremove,:) = [];
tp@0 625 expandedivsd(ivremove) = [];
tp@0 626 nposs = length(expandedivsd);
tp@0 627
tp@0 628 test = [expandedposscombs expandedpossedges];
tp@0 629 ncombs = length(expandedpossedges);
tp@0 630 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 631 ivremove = find(dtest==1);
tp@0 632
tp@0 633 end
tp@0 634
tp@0 635 if SHOWTEXT >= 3
tp@0 636 disp([' ',int2str(nposs),' IS+edge combinations after edge segments are packed together'])
tp@0 637 end
tp@0 638 combstoremove = setdiff(ivsd,expandedivsd);
tp@0 639 addtocomb(combstoremove,:) = [];
tp@0 640 addtooriginsfrom(combstoremove) = [];
tp@0 641 addtoISESVISIBILITY(combstoremove) = [];
tp@0 642 nadds = length(addtooriginsfrom);
tp@0 643 end
tp@0 644
tp@0 645 end
tp@0 646
tp@0 647 % Combinations of diffraction-diffraction don't need to be checked because
tp@0 648 % 'edgeseesedge' should have taken care of this.
tp@0 649
tp@0 650 %----------------------------------------------------------------------
tp@0 651 % Add the new combinations to the master list
tp@0 652
tp@0 653 if nadds > 0
tp@0 654 POTENTIALISES = [POTENTIALISES;[addtocomb zeros(nadds,specorder-2)]];
tp@0 655 ORIGINSFROM = [ORIGINSFROM;addtooriginsfrom];
tp@0 656
tp@0 657 if difforder > 0
tp@0 658 ivtemp = find(addtocomb(:,1)<=nplanes & addtocomb(:,2)>nplanes);
tp@0 659 if ~isempty(ivtemp),
tp@0 660 addtoISESVISIBILITY(ivtemp) = edgeweightlist;
tp@0 661 end
tp@0 662 ivtemp = find(addtocomb(:,1)>nplanes);
tp@0 663 end
tp@0 664 ISESVISIBILITY = [ISESVISIBILITY;addtoISESVISIBILITY];
tp@0 665
tp@0 666 endindices(2) = length(ORIGINSFROM);
tp@0 667
tp@0 668 % Compute the IS coords of the combinations with only specular reflections
tp@0 669
tp@0 670 ISCOORDStoadd = zeros(nadds,3);
tp@0 671 if difforder > 0
tp@0 672 ivss = uint32(find(addtocomb(:,1)<=nplanes & addtocomb(:,2)<=nplanes));
tp@0 673 else
tp@0 674 ivss = uint32([1:nadds].');
tp@0 675 end
tp@0 676 soutoreflect = ISCOORDS(ORIGINSFROM(double(ivss)+endindices(1)),1:3);
tp@0 677 ISCOORDStoadd(ivss,:) = EDB1findis(soutoreflect,POTENTIALISES(double(ivss)+endindices(1),2),planeeqs,size(soutoreflect,1),onesvec3);
tp@0 678 clear soutoreflect
tp@0 679 ISCOORDS = [ISCOORDS;ISCOORDStoadd];
tp@0 680
tp@0 681 end
tp@0 682 end
tp@0 683 end
tp@0 684 end
tp@0 685
tp@0 686 %##################################################################
tp@0 687 %##################################################################
tp@0 688 %##################################################################
tp@0 689 %
tp@0 690 % Third order and higher
tp@0 691 %
tp@0 692 %---------------------------------------------------------------------------------------------
tp@0 693
tp@0 694 for ordernum = 3:specorder
tp@0 695 if SHOWTEXT >= 3
tp@0 696 disp([' Order number ',int2str(ordernum)])
tp@0 697 end
tp@0 698
tp@0 699 startindices(ordernum) = startindices(ordernum-1) + length(addtocomb);
tp@0 700
tp@0 701 % The lines below could use
tp@0 702 % planestoprop = unique(addtocomb(:,ordernum-1));
tp@0 703
tp@0 704 planelist = sort(addtocomb(:,ordernum-1));
tp@0 705
tp@0 706 if ~isempty(planelist),
tp@0 707 planestoprop = planelist([1;find(diff(double(planelist))~=0)+1]);
tp@0 708 else
tp@0 709 planestoprop = [];
tp@0 710 end
tp@0 711 clear planelist
tp@0 712
tp@0 713 maxnumberofvispla = max(sum(PLANEseesPLANE(:,planestoprop)==1));
tp@0 714 oldaddtocomb = addtocomb;
tp@0 715
tp@0 716 nadds = size(addtocomb,1);
tp@0 717 if nhighval < 255
tp@0 718 addtocomb = uint8(zeros(nadds*maxnumberofvispla,ordernum));
tp@0 719 else
tp@0 720 addtocomb = uint16(zeros(nadds*maxnumberofvispla,ordernum));
tp@0 721 end
tp@0 722
tp@0 723 onesvec2 = ones(1,maxnumberofvispla);
tp@0 724 for ii = 1:ordernum-1
tp@0 725 temp = oldaddtocomb(:,ii);
tp@0 726 addtocomb(:,ii) = reshape(temp(:,onesvec2).',length(temp)*maxnumberofvispla,1);
tp@0 727 end
tp@0 728 nrows = size(addtocomb,1);
tp@0 729 clear oldaddtocomb temp
tp@0 730
tp@0 731 if nrows > 0
tp@0 732 startpos = [[1:maxnumberofvispla:nrows] nrows+1 ].';
tp@0 733 else
tp@0 734 startpos = 1;
tp@0 735 end
tp@0 736
tp@0 737 for ii = 1:length(startpos)-1
tp@0 738 if nvisPLANES(addtocomb(startpos(ii),ordernum-1)) > 0
tp@0 739 possibleplanes = listofvisPLANES(addtocomb(startpos(ii),ordernum-1),1:nvisPLANES(addtocomb(startpos(ii),ordernum-1))).';
tp@0 740 nposs = length(possibleplanes);
tp@0 741 addtocomb( startpos(ii):startpos(ii)+nposs-1,ordernum) = possibleplanes;
tp@0 742 end
tp@0 743 end
tp@0 744
tp@0 745 addtooriginsfrom = uint32([1:nadds].' + startindices(ordernum-1)-1);
tp@0 746 addtooriginsfrom = addtooriginsfrom(:,ones(1,maxnumberofvispla));
tp@0 747 addtooriginsfrom = reshape(addtooriginsfrom.',maxnumberofvispla*nadds,1);
tp@0 748
tp@0 749 addtoISESVISIBILITY = reshape(addtoISESVISIBILITY(:,ones(1,maxnumberofvispla)).',nadds*maxnumberofvispla,1);
tp@0 750
tp@0 751 clear startpos
tp@0 752 cutout = find(addtocomb(:,ordernum)==0);
tp@0 753 addtocomb(cutout,:) = [];
tp@0 754 addtooriginsfrom(cutout) = [];
tp@0 755 addtoISESVISIBILITY(cutout) = [];
tp@0 756 nadds = size(addtocomb,1);
tp@0 757
tp@0 758 if ordernum == specorder
tp@0 759 clear cutout
tp@0 760 end
tp@0 761
tp@0 762 %--------------------------------------------------------------------------
tp@0 763 % Try to get rid of some combinations that we know can not be propagated
tp@0 764
tp@0 765 % Find those combinations of specular-specular where the IS is behind the second reflecting plane
tp@0 766 % because they can then never ever be propagated.
tp@0 767
tp@0 768 if difforder > 0
tp@0 769 ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0 770 else
tp@0 771 ivss = uint32([1:nadds].');
tp@0 772 end
tp@0 773
tp@0 774 imsoucheck = dot((ISCOORDS(addtooriginsfrom(ivss),1:3) - corners(planecorners(addtocomb(ivss,ordernum),1),:)).',planenvecs(addtocomb(ivss,ordernum),:).').';
tp@0 775 GEOMACC = 1e-10;
tp@0 776 imsounotOK = uint32(find(imsoucheck < GEOMACC));
tp@0 777 clear imsoucheck
tp@0 778 addtocomb(ivss(imsounotOK),:) = [];
tp@0 779 addtooriginsfrom(ivss(imsounotOK),:) = [];
tp@0 780 addtoISESVISIBILITY(ivss(imsounotOK)) = [];
tp@0 781 clear imsounotOK
tp@0 782
tp@0 783 % Addition 20100210 PS: we should check all combinations with dss in the
tp@0 784 % three last orders of addtocomb: if the edge in the first 'd' belongs to the
tp@0 785 % plane of the last 's' *and* the two 's' planes form a 90 degree
tp@0 786 % corner, then these combos should be removed.
tp@0 787 %
tp@0 788 % As a first step, we check if there are any interior-90-degree
tp@0 789 % corners/edges of the model, because only then can these combos occur.
tp@0 790
tp@0 791 if ordernum == 3
tp@0 792 deg90edges = find( abs(closwedangvec-3*pi/2) < GEOMACC );
tp@0 793 if any(deg90edges)
tp@0 794 deg90planepairs = planesatedge(deg90edges,:);
tp@0 795 end
tp@0 796 end
tp@0 797
tp@0 798 if difforder > 0 & any(deg90edges)
tp@0 799 ivdss = uint32(find( double(addtocomb(:,ordernum-2)>nplanes).*double(addtocomb(:,ordernum-1)<=nplanes).*double(addtocomb(:,ordernum)<=nplanes) ));
tp@0 800 else
tp@0 801 ivdss = [];
tp@0 802 end
tp@0 803
tp@0 804 if ~isempty(ivdss)
tp@0 805
tp@0 806 % First we check if the 'd' (the edge) belongs to the last 's'
tp@0 807 A = addtocomb(ivdss,ordernum-2:ordernum);
tp@0 808
tp@0 809 ivsubset = find( (planesatedge(A(:,1)-nplanes,1) == A(:,3)) | (planesatedge(A(:,1)-nplanes,2) == A(:,3)) );
tp@0 810 ivdss = ivdss(ivsubset);
tp@0 811
tp@0 812 if ~isempty(ivdss)
tp@0 813 % Second, we check if the two 's' planes form a 90 deg corner
tp@0 814 planesare90degpairs = ismember( A(:,2),deg90planepairs(:,1) ).*ismember( A(:,3),deg90planepairs(:,2) ) + ismember( A(:,2),deg90planepairs(:,2) ).*ismember( A(:,3),deg90planepairs(:,1) );
tp@0 815 ivsubset = find(planesare90degpairs);
tp@0 816 ivdss = ivdss(ivsubset);
tp@0 817
tp@0 818 if ~isempty(ivdss)
tp@0 819 addtocomb(ivdss,:) = [];
tp@0 820 addtooriginsfrom(ivdss,:) = [];
tp@0 821 addtoISESVISIBILITY(ivdss) = [];
tp@0 822 clear ivdss ivsubset planesare90degpairs
tp@0 823 end
tp@0 824 end
tp@0 825 end
tp@0 826
tp@0 827 nadds = size(addtocomb,1);
tp@0 828
tp@0 829 if nadds > 0 & doconetrace == 1
tp@0 830
tp@0 831 if difforder > 0
tp@0 832 ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0 833 else
tp@0 834 ivss = uint32([1:nadds].');
tp@0 835 end
tp@0 836
tp@0 837 beamdirection1 = planemidpoints(addtocomb(ivss,ordernum-1),:) - ISCOORDS(addtooriginsfrom(ivss),1:3);
tp@0 838 beamlength = sqrt( sum( beamdirection1.'.^2 ) ).';
tp@0 839 beamdirection1 = beamdirection1./beamlength(:,ones(1,3));
tp@0 840 maxprojradius = maxdisttocorner(addtocomb(ivss,ordernum-1));
tp@0 841 iv = find(beamlength<maxprojradius);
tp@0 842 if ~isempty(iv)
tp@0 843 beamlength(iv) = beamlength(iv)*0+ eps*1e6;
tp@0 844 end
tp@0 845 beamangle1 = atan(maxprojradius./beamlength);
tp@0 846
tp@0 847 % Now we calculate the beam-widths to the second reflection
tp@0 848 % plane, seen from the first-order image source. If this
tp@0 849 % beamwidth is outside the previously calculated beam of
tp@0 850 % reflection plane 1, then plane 2 can never be seen
tp@0 851 % through plane 1.
tp@0 852
tp@0 853 beamdirection2 = planemidpoints(addtocomb(ivss,ordernum),:) - ISCOORDS(addtooriginsfrom(ivss),1:3);
tp@0 854 beamlength = sqrt( sum( beamdirection2.'.^2 ) ).';
tp@0 855 beamdirection2 = beamdirection2./beamlength(:,ones(1,3));
tp@0 856 maxprojradius = maxdisttocorner(addtocomb(ivss,ordernum));
tp@0 857 iv = find(beamlength<maxprojradius);
tp@0 858 if ~isempty(iv)
tp@0 859 beamlength(iv) = beamlength(iv)*0+ eps*10;
tp@0 860 end
tp@0 861 beamangle2 = atan(maxprojradius./beamlength);
tp@0 862 clear maxprojradius beamlength
tp@0 863
tp@0 864 ivcutout = find( (beamangle1 + beamangle2 < acos(sum( (beamdirection1.*beamdirection2).' ).')) & (beamangle1 <1.56) & (beamangle2 < 1.56));
tp@0 865 clear beamdirection1 beamdirection2 beamangle1 beamangle2
tp@0 866
tp@0 867 if ~isempty(ivcutout)
tp@0 868 addtocomb(ivss(ivcutout),:) = [];
tp@0 869 addtooriginsfrom(ivss(ivcutout),:) = [];
tp@0 870 addtoISESVISIBILITY(ivss(ivcutout)) = [];
tp@0 871 nadds = size(addtocomb,1);
tp@0 872 end
tp@0 873
tp@0 874 end
tp@0 875
tp@0 876 % Combinations with all-specular plus diffraction as the last one:
tp@0 877 % check if the IS is behind both of the two planes that the reflecting edge connects to.
tp@0 878 %
tp@0 879 % Bug found 080711: The visibility test must be split into two, depending on the wedge angle!!
tp@0 880 % If the open wedge angle > pi, then it is correct to check if the IS behind both of the planes.
tp@0 881 % If the open wedge angle < pi, then it should be checked if the IS behind one of the planes!!
tp@0 882
tp@0 883 if difforder > 0
tp@0 884
tp@0 885 ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0 886
tp@0 887 if ~isempty(ivsd)
tp@0 888
tp@0 889 % First we check the sssssd combinations where the open wedge angle > pi
tp@0 890 % For those cases, we can exclude combinations where the IS
tp@0 891 % is behind both planes
tp@0 892
tp@0 893 ivsd1 = ivsd( find( closwedangvec( double(addtocomb(ivsd,ordernum))-nplanes ) < pi ));
tp@0 894
tp@0 895 if ~isempty(ivsd1)
tp@0 896 imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1) ,:).').';
tp@0 897 imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2) ,:).').';
tp@0 898 GEOMACC = 1e-10;
tp@0 899 imsounotOK = find(imsoucheck1 < GEOMACC & imsoucheck2 < GEOMACC);
tp@0 900 clear imsoucheck1 imsoucheck2
tp@0 901 addtocomb(ivsd1(imsounotOK),:) = [];
tp@0 902 addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0 903 addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0 904 clear imsounotOK
tp@0 905 end
tp@0 906
tp@0 907 % Second we check the sssssd combinations where the open wedge angle < pi
tp@0 908 % For those cases, we can exclude combinations where the IS
tp@0 909 % is behind one of the two planes
tp@0 910
tp@0 911 ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0 912 if ~isempty(ivsd)
tp@0 913 ivsd1 = ivsd( find( closwedangvec( double(addtocomb(ivsd,ordernum))-nplanes ) > pi ));
tp@0 914
tp@0 915 if ~isempty(ivsd1)
tp@0 916 imsoucheck1 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,1) ,:).').';
tp@0 917 imsoucheck2 = dot((ISCOORDS(addtooriginsfrom(ivsd1),1:3) - corners(planecorners( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2) ,1),:)).',planenvecs( planesatedge(double(addtocomb(ivsd1,ordernum))-nplanes,2) ,:).').';
tp@0 918 GEOMACC = 1e-10;
tp@0 919 imsounotOK = find(imsoucheck1 < GEOMACC | imsoucheck2 < GEOMACC);
tp@0 920 clear imsoucheck1 imsoucheck2
tp@0 921 addtocomb(ivsd1(imsounotOK),:) = [];
tp@0 922 addtooriginsfrom(ivsd1(imsounotOK),:) = [];
tp@0 923 addtoISESVISIBILITY(ivsd1(imsounotOK)) = [];
tp@0 924 clear imsounotOK
tp@0 925 end
tp@0 926 end
tp@0 927
tp@0 928 % For combinations including dsd, if the two diff. edges are
tp@0 929 % aligned and the intermediate specular reflections are
tp@0 930 % caused by planes that are perpendicular to the edges
tp@0 931 % remove these combinations
tp@0 932
tp@0 933 if difforder >= 2
tp@0 934 for kk = 1:ordernum-2
tp@0 935 ncombs = size(addtocomb,1);
tp@0 936 matchvec = ones(ncombs,1);
tp@0 937 nprespecs = kk-1;
tp@0 938 nmidspecs = ordernum-1-kk;
tp@0 939 matchvec = matchvec.*(addtocomb(:,ordernum)>nplanes).*(addtocomb(:,kk)>nplanes);
tp@0 940 if nmidspecs == 1
tp@0 941 matchvec = matchvec.*(addtocomb(:,2)<=nplanes);
tp@0 942 elseif nmidspecs >= 2
tp@0 943 matchvec = matchvec.*prod( double(addtocomb(:,kk+1:ordernum-1)<=nplanes).' ).';
tp@0 944 end
tp@0 945 ivdsd = uint32(find( matchvec ));
tp@0 946 if ~isempty(ivdsd)
tp@0 947 edge1 = double(addtocomb(ivdsd,kk)) - nplanes;
tp@0 948 edge2 = double(addtocomb(ivdsd,ordernum)) - nplanes;
tp@0 949 midspecs = double(addtocomb(ivdsd,kk+1:ordernum-1));
tp@0 950
tp@0 951 lookupindmat1 = (edge1-1)*nedges + edge2;
tp@0 952 matrixcolnumbers = (edge1-1)*nplanes;
tp@0 953 lookupindmat2 = matrixcolnumbers(:,ones(1,nmidspecs)) + midspecs;
tp@0 954
tp@0 955 if nmidspecs == 1
tp@0 956 ivinvalid = find(edgealignedwithedge(lookupindmat1) & edgeperptoplane(lookupindmat2));
tp@0 957 else
tp@0 958 ivinvalid = find(edgealignedwithedge(lookupindmat1) & prod(edgeperptoplane(lookupindmat2).').');
tp@0 959 end
tp@0 960 addtocomb(ivdsd(ivinvalid),:) = [];
tp@0 961 addtooriginsfrom(ivdsd(ivinvalid),:) = [];
tp@0 962 addtoISESVISIBILITY(ivdsd(ivinvalid)) = [];
tp@0 963 end
tp@0 964 end
tp@0 965 end
tp@0 966
tp@0 967 % We need to find the valid combinations again after we have removed a
tp@0 968 % number of combinations
tp@0 969
tp@0 970 ivsd = uint32(find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes)));
tp@0 971
tp@0 972 %----------------------------------------------------------------------
tp@0 973 % Combinations of specular-diffraction that are not visible at
tp@0 974 % all should be removed and then they are not propagated either.
tp@0 975
tp@0 976 if nhighval < 255
tp@0 977 possibleedges = uint8(double(addtocomb(ivsd,ordernum)) - nplanes);
tp@0 978 possiblecombs = uint8(double(addtocomb(ivsd,1:ordernum-1)));
tp@0 979 else
tp@0 980 possibleedges = uint16(double(addtocomb(ivsd,ordernum)) - nplanes);
tp@0 981 possiblecombs = uint16(double(addtocomb(ivsd,1:ordernum-1)));
tp@0 982 end
tp@0 983 reftoISCOORDS = addtooriginsfrom(ivsd);
tp@0 984
tp@0 985 % Expand to take the edge subdivisions into account
tp@0 986
tp@0 987 nposs = length(ivsd);
tp@0 988 nposs = nposs*nedgesubs; % We treat all the little edge subdivisions as separate edges
tp@0 989
tp@0 990 expandedposscombs = reshape(repmat(possiblecombs.',[nedgesubs,1]),ordernum-1,nposs).';
tp@0 991 clear possiblecombs
tp@0 992 expandedreftoISCOORDS = reftoISCOORDS(:,onesvec1);
tp@0 993 expandedreftoISCOORDS = reshape(expandedreftoISCOORDS.',nposs,1);
tp@0 994 expandedpossedges = possibleedges(:,onesvec1);
tp@0 995 expandedpossedges = reshape(expandedpossedges.',nposs,1);
tp@0 996 expandedivsd = ivsd(:,onesvec1);
tp@0 997 expandedivsd = reshape(expandedivsd.',nposs,1);
tp@0 998
tp@0 999 if SHOWTEXT >= 3
tp@0 1000 disp([' ',int2str(nposs),' IS+edge segments found initially '])
tp@0 1001 end
tp@0 1002
tp@0 1003 % Go through, iteratively, and check if the path from S to edge passes
tp@0 1004 % through all reflection planes along the way
tp@0 1005
tp@0 1006 for jj = ordernum-1:-1:1
tp@0 1007
tp@0 1008 if nposs > 0
tp@0 1009
tp@0 1010 if jj == ordernum-1
tp@0 1011 fromcoords = ISCOORDS(expandedreftoISCOORDS,:);
tp@0 1012 [tocoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0 1013 clear possibleedges
tp@0 1014 else
tp@0 1015 eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])
tp@0 1016 ivlist = ORIGINSFROM(expandedreftoISCOORDS);
tp@0 1017 for kk = jj:ordernum-3
tp@0 1018 ivlist = ORIGINSFROM(ivlist);
tp@0 1019 end
tp@0 1020 fromcoords = ISCOORDS(ivlist,:);
tp@0 1021
tp@0 1022 end
tp@0 1023
tp@0 1024 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,jj),4),planenvecs(expandedposscombs(:,jj),:),minvals(expandedposscombs(:,jj),:),...
tp@0 1025 maxvals(expandedposscombs(:,jj),:),planecorners(expandedposscombs(:,jj),:),corners,ncornersperplanevec(expandedposscombs(:,jj)));
tp@0 1026 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 1027 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 1028 disp(' handled correctly yet.')
tp@0 1029 end
tp@0 1030 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0 1031
tp@0 1032 expandedivsd = expandedivsd(hitplanes);
tp@0 1033 expandedposscombs = expandedposscombs(hitplanes,:);
tp@0 1034 expandedreftoISCOORDS = expandedreftoISCOORDS(hitplanes);
tp@0 1035 expandedpossedges = expandedpossedges(hitplanes);
tp@0 1036 edgeweightlist = edgeweightlist(hitplanes);
tp@0 1037 toedgecoords = tocoords(hitplanes,:);
tp@0 1038 if jj < ordernum-1
tp@0 1039 for kk = jj+1:ordernum-1
tp@0 1040 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);'])
tp@0 1041 end
tp@0 1042 end
tp@0 1043
tp@0 1044 nposs = length(expandedivsd);
tp@0 1045
tp@0 1046 end
tp@0 1047 if SHOWTEXT >= 3
tp@0 1048 disp([' ',int2str(nposs),' IS+edge segments survived the visibility test in refl plane ',JJ(jj,1:JJnumbofchars(jj))])
tp@0 1049 end
tp@0 1050
tp@0 1051 end
tp@0 1052
tp@0 1053 if obstructtestneeded & nposs > 0
tp@0 1054 for jj = 1:ordernum
tp@0 1055 if nposs > 0
tp@0 1056
tp@0 1057 if jj==1
tp@0 1058 fromcoords = S;
tp@0 1059 startplanes = [];
tp@0 1060 else
tp@0 1061 startplanes = expandedposscombs(:,jj-1);
tp@0 1062 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0 1063 end
tp@0 1064 if jj == ordernum
tp@0 1065 tocoords = toedgecoords;
tp@0 1066 endplanes = [planesatedge(expandedpossedges,1) planesatedge(expandedpossedges,2)];
tp@0 1067 else
tp@0 1068 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
tp@0 1069 endplanes = expandedposscombs(:,jj);
tp@0 1070 end
tp@0 1071
tp@0 1072 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 1073 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 1074
tp@0 1075 if nobstructions > 0
tp@0 1076 expandedivsd = expandedivsd(nonobstructedpaths);
tp@0 1077 expandedposscombs = expandedposscombs(nonobstructedpaths,:);
tp@0 1078 expandedreftoISCOORDS = expandedreftoISCOORDS(nonobstructedpaths);
tp@0 1079 expandedpossedges = expandedpossedges(nonobstructedpaths);
tp@0 1080 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 1081 toedgecoords = tocoords(nonobstructedpaths,:);
tp@0 1082 nposs = length(expandedivsd);
tp@0 1083 for kk = 1:ordernum-1
tp@0 1084 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])
tp@0 1085 end
tp@0 1086
tp@0 1087 end
tp@0 1088 end
tp@0 1089 end
tp@0 1090 end
tp@0 1091
tp@0 1092 if SHOWTEXT >= 3
tp@0 1093 disp([' ',int2str(nposs),' IS+edge segments survived the obstruction test'])
tp@0 1094 end
tp@0 1095
tp@0 1096 % There are repetitions of each combination since each edge segment has
tp@0 1097 % been treated separately. Pack them together now
tp@0 1098
tp@0 1099 test = [expandedposscombs expandedpossedges];
tp@0 1100 ncombs = length(expandedpossedges);
tp@0 1101 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 1102 ivremove = find(dtest==1);
tp@0 1103
tp@0 1104 while ~isempty(ivremove)
tp@0 1105
tp@0 1106 edgeweightlist(ivremove+1) = double(edgeweightlist(ivremove+1)) + double(edgeweightlist(ivremove));
tp@0 1107 edgeweightlist(ivremove) = [];
tp@0 1108 expandedpossedges(ivremove) = [];
tp@0 1109 expandedposscombs(ivremove,:) = [];
tp@0 1110 expandedivsd(ivremove) = [];
tp@0 1111 nposs = length(expandedivsd);
tp@0 1112
tp@0 1113 test = [expandedposscombs expandedpossedges];
tp@0 1114 ncombs = length(expandedpossedges);
tp@0 1115 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 1116 ivremove = find(dtest==1);
tp@0 1117
tp@0 1118 end
tp@0 1119
tp@0 1120 if SHOWTEXT >= 3
tp@0 1121 disp([' ',int2str(nposs),' IS+edge segments after edge segments are packed together'])
tp@0 1122 end
tp@0 1123
tp@0 1124 combstoremove = setdiff(ivsd,expandedivsd);
tp@0 1125 addtocomb(combstoremove,:) = [];
tp@0 1126 addtooriginsfrom(combstoremove) = [];
tp@0 1127 addtoISESVISIBILITY(combstoremove) = [];
tp@0 1128 nadds = length(addtooriginsfrom);
tp@0 1129 end
tp@0 1130 end
tp@0 1131
tp@0 1132 POTENTIALISES = [POTENTIALISES;[addtocomb zeros(nadds,specorder-ordernum)]];
tp@0 1133 ORIGINSFROM = [ORIGINSFROM;addtooriginsfrom];
tp@0 1134
tp@0 1135 ivtemp = find(prod( double(addtocomb(:,1:ordernum-1)<=nplanes).' ).' & (addtocomb(:,ordernum)>nplanes));
tp@0 1136 if difforder > 0
tp@0 1137 if ~isempty(ivtemp)
tp@0 1138 addtoISESVISIBILITY(ivtemp) = edgeweightlist;
tp@0 1139 end
tp@0 1140 end
tp@0 1141 ISESVISIBILITY = [ISESVISIBILITY;addtoISESVISIBILITY];
tp@0 1142
tp@0 1143 endindices(ordernum) = length(ORIGINSFROM);
tp@0 1144
tp@0 1145 % Compute the IS coords
tp@0 1146
tp@0 1147 ISCOORDStoadd = zeros(nadds,3);
tp@0 1148 if difforder > 0
tp@0 1149 ivss = uint32(find(prod( double(addtocomb<=nplanes).' ).'));
tp@0 1150 else
tp@0 1151 ivss = uint32([1:nadds].');
tp@0 1152 end
tp@0 1153 soutoreflect = ISCOORDS(ORIGINSFROM(double(ivss)+endindices(ordernum-1)),1:3);
tp@0 1154 ISCOORDStoadd(ivss,:) = EDB1findis(soutoreflect,POTENTIALISES(double(ivss)+endindices(ordernum-1),ordernum),planeeqs,size(soutoreflect,1),onesvec3);
tp@0 1155 clear soutoreflect
tp@0 1156 ISCOORDS = [ISCOORDS;ISCOORDStoadd];
tp@0 1157 clear ISCOORDStoadd
tp@0 1158
tp@0 1159 end
tp@0 1160
tp@0 1161 ntot = size(POTENTIALISES,1);
tp@0 1162
tp@0 1163 %-------------------------------------------------------
tp@0 1164 % Add the direct sound to the start of the list
tp@0 1165
tp@0 1166 POTENTIALISES = [zeros(1,specorder);POTENTIALISES];
tp@0 1167 ORIGINSFROM = [0;ORIGINSFROM];
tp@0 1168 startindices = startindices+1;
tp@0 1169 endindices = endindices+1;
tp@0 1170 ORIGINSFROM = uint32(double(ORIGINSFROM)+1);
tp@0 1171
tp@0 1172 ISCOORDS = [S;ISCOORDS];
tp@0 1173 ISESVISIBILITY = [0;ISESVISIBILITY];
tp@0 1174
tp@0 1175 %-------------------------------------------------------
tp@0 1176 % Trim the list if there are zeros-only in the last column(s)
tp@0 1177
tp@0 1178 [n1,n2] = size(POTENTIALISES);
tp@0 1179 if n2 > 1
tp@0 1180 columnstokeep = find(sum(POTENTIALISES)~=0);
tp@0 1181 POTENTIALISES = POTENTIALISES(:,columnstokeep);
tp@0 1182 [n1,n2new] = size(POTENTIALISES);
tp@0 1183 if n2new < n2
tp@0 1184 specorderold = specorder;
tp@0 1185 specorder = n2new;
tp@0 1186 end
tp@0 1187 end
tp@0 1188
tp@0 1189 %-------------------------------------------------------
tp@0 1190 % Create index lists
tp@0 1191
tp@0 1192 % First, the total reflection order
tp@0 1193
tp@0 1194 [n1,n2] = size(POTENTIALISES);
tp@0 1195
tp@0 1196 lengthNspecmatrix = [];
tp@0 1197 lengthNdiffmatrix = [];
tp@0 1198 singlediffcol = [];
tp@0 1199 startindicessinglediff = [];
tp@0 1200 endindicessinglediff = [];
tp@0 1201
tp@0 1202 if n1 > 0 & n2 > 0
tp@0 1203 if n2 > 1
tp@0 1204 REFLORDER = uint8(sum(POTENTIALISES.'>0).');
tp@0 1205 else
tp@0 1206 REFLORDER = uint8(POTENTIALISES>0);
tp@0 1207 end
tp@0 1208
tp@0 1209 if difforder > 0
tp@0 1210 if specorder > 1
tp@0 1211 ivss = uint32(find(prod( double(POTENTIALISES<=nplanes).' ).'));
tp@0 1212 else
tp@0 1213 ivss = uint32(find( POTENTIALISES<=nplanes ));
tp@0 1214 end
tp@0 1215
tp@0 1216 ivother = uint32([1:length(ORIGINSFROM)].');
tp@0 1217 ivother(ivss) = [];
tp@0 1218 ivss(1) = [];
tp@0 1219 IVNDIFFMATRIX = uint32(zeros(2,specorder));
tp@0 1220 lengthNdiffmatrix = zeros(1,specorder);
tp@0 1221 nrowsdiff = 2;
tp@0 1222 IVNSPECMATRIX = uint32(zeros(2,specorder));
tp@0 1223 lengthNspecmatrix = zeros(1,specorder);
tp@0 1224 nrowsspec = 2;
tp@0 1225
tp@0 1226 for ii = 1:specorder
tp@0 1227 if specorder > 1
tp@0 1228 ivreftoNdiff = uint32(find(sum( double(POTENTIALISES(ivother,:)> nplanes).' ).'==ii));
tp@0 1229 else
tp@0 1230 ivreftoNdiff = uint32(find(( double(POTENTIALISES(ivother,:)> nplanes).' ).'==ii));
tp@0 1231 end
tp@0 1232 ivreftoNspec = uint32(find(REFLORDER(ivss)==ii));
tp@0 1233
tp@0 1234 ivNdiff = ivother(ivreftoNdiff);
tp@0 1235 if ~isempty(ivNdiff)
tp@0 1236 lengthNdiffmatrix(ii) = length(ivNdiff);
tp@0 1237 if lengthNdiffmatrix(ii) > nrowsdiff
tp@0 1238 IVNDIFFMATRIX = [IVNDIFFMATRIX;uint32(zeros(lengthNdiffmatrix(ii)-nrowsdiff,specorder))];
tp@0 1239 nrowsdiff = lengthNdiffmatrix(ii);
tp@0 1240 end
tp@0 1241 IVNDIFFMATRIX(1: lengthNdiffmatrix(ii),ii) = ivNdiff;
tp@0 1242 end
tp@0 1243 ivother(ivreftoNdiff) = [];
tp@0 1244
tp@0 1245 ivNspec = ivss(ivreftoNspec);
tp@0 1246 if ~isempty(ivNspec)
tp@0 1247 lengthNspecmatrix(ii) = length(ivNspec);
tp@0 1248 if lengthNspecmatrix(ii) > nrowsspec
tp@0 1249 IVNSPECMATRIX = [IVNSPECMATRIX;uint32(zeros(lengthNspecmatrix(ii)-nrowsspec,specorder))];
tp@0 1250 nrowsspec = lengthNspecmatrix(ii);
tp@0 1251 end
tp@0 1252 IVNSPECMATRIX(1: lengthNspecmatrix(ii),ii) = ivNspec;
tp@0 1253 end
tp@0 1254 end
tp@0 1255
tp@0 1256 % Determine which column the single diffraction is in
tp@0 1257
tp@0 1258 test = POTENTIALISES(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1),:)>nplanes;
tp@0 1259
tp@0 1260 colweight = [1:specorder];
tp@0 1261 nsingles = length(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1));
tp@0 1262 if specorder > 1
tp@0 1263 singlediffcol = uint8(sum( (test.*colweight(ones(nsingles,1),:)).').');
tp@0 1264 else
tp@0 1265 singlediffcol = uint8(ones(nsingles,1));
tp@0 1266 end
tp@0 1267
tp@0 1268 startindicessinglediff = zeros(specorder,1);
tp@0 1269 endindicessinglediff = zeros(specorder,1);
tp@0 1270 startindicessinglediff(1) = 1;
tp@0 1271 for ii = 1:specorder-1
tp@0 1272 iv = find(IVNDIFFMATRIX(1:lengthNdiffmatrix(1),1)<=endindices(ii));
tp@0 1273 if ~isempty(iv)
tp@0 1274 endindicessinglediff(ii) = iv(length(iv));
tp@0 1275 else
tp@0 1276 endindicessinglediff(ii) = 0;
tp@0 1277 end
tp@0 1278 startindicessinglediff(ii+1) = endindicessinglediff(ii)+1;
tp@0 1279 end
tp@0 1280 endindicessinglediff(specorder) = lengthNdiffmatrix(1);
tp@0 1281 else
tp@0 1282 ivss = uint32([2:endindices(length(endindices))].');
tp@0 1283 IVNSPECMATRIX = uint32(zeros(2,specorder));
tp@0 1284 lengthNspecmatrix = zeros(1,specorder);
tp@0 1285 nrowsspec = 2;
tp@0 1286
tp@0 1287 for ii = 1:specorder
tp@0 1288 ivreftoNspec = uint32(find(REFLORDER(ivss)==ii));
tp@0 1289
tp@0 1290 ivNspec = ivss(ivreftoNspec);
tp@0 1291 lengthNspecmatrix(ii) = length(ivNspec);
tp@0 1292 if lengthNspecmatrix(ii) > nrowsspec
tp@0 1293 IVNSPECMATRIX = [IVNSPECMATRIX;uint32(zeros(lengthNspecmatrix(ii)-nrowsspec,specorder))];
tp@0 1294 nrowsspec = lengthNspecmatrix(ii);
tp@0 1295 end
tp@0 1296 IVNSPECMATRIX(1: lengthNspecmatrix(ii),ii) = ivNspec;
tp@0 1297
tp@0 1298 end
tp@0 1299
tp@0 1300 IVNDIFFMATRIX = [];
tp@0 1301 lengthNdiffmatrix = [];
tp@0 1302
tp@0 1303 singlediffcol = [];
tp@0 1304 startindicessinglediff = [];
tp@0 1305 endindicessinglediff = [];
tp@0 1306 end
tp@0 1307
tp@0 1308 end
tp@0 1309
tp@0 1310 %-------------------------------------------------------
tp@0 1311 % Create a pointer list so that it is possible to find a combination
tp@0 1312 % that ends with edge-plane2, given an edge-plane1-plane2 combination
tp@0 1313 %
tp@0 1314 % NB!! The list points to the original POTENTIALISES (and related lists)
tp@0 1315
tp@0 1316 if difforder > 0 & specorder > 1
tp@0 1317
tp@0 1318 if SHOWTEXT >= 3
tp@0 1319 disp([' Building a pointer list for image receiver combinations'])
tp@0 1320 end
tp@0 1321
tp@0 1322 % First select only those that have a single diff combo and that
tp@0 1323 % start with the diffraction
tp@0 1324
tp@0 1325 ndecimaldivider = (nplanes+2);
tp@0 1326 PointertoIRcombs = sparse(zeros(ndecimaldivider^(specorder-1),1));
tp@0 1327 IRoriginsfrom = uint32(zeros(size(ORIGINSFROM)));
tp@0 1328
tp@0 1329 for ii = 1:specorder-1
tp@0 1330 if SHOWTEXT >= 3
tp@0 1331 disp([' order ',int2str(ii)])
tp@0 1332 end
tp@0 1333
tp@0 1334 if ii == 1
tp@0 1335 ivrange = uint32([startindicessinglediff(2):endindicessinglediff(2)]);
tp@0 1336 masterivlistselect = IVNDIFFMATRIX(ivrange,1);
tp@0 1337 masterivlistselect = masterivlistselect(find(singlediffcol(ivrange)==1));
tp@0 1338
tp@0 1339 % IRoriginsfrom should be given the value zero for these
tp@0 1340 % first-order specular combos (after one diff) so we don't
tp@0 1341 % bother assigning any value
tp@0 1342
tp@0 1343 [B,IA,JA] = unique(POTENTIALISES(masterivlistselect,2));
tp@0 1344 PointertoIRcombs(B) = masterivlistselect(IA);
tp@0 1345
tp@0 1346 % Now we check if there any active wall numbers that don't
tp@0 1347 % occur in an edge-spec combination. For those combinations
tp@0 1348 % we can point to a specular combos instead, that ends with the
tp@0 1349 % same plane number.
tp@0 1350
tp@0 1351 wallist = POTENTIALISES(IVNSPECMATRIX(1:lengthNspecmatrix(1),1),1);
tp@0 1352 ivnotincludedyet = find(~ismember(wallist,B));
tp@0 1353 if ~isempty(ivnotincludedyet)
tp@0 1354 PointertoIRcombs(wallist(ivnotincludedyet)) = IVNSPECMATRIX(ivnotincludedyet,1);
tp@0 1355 end
tp@0 1356
tp@0 1357 elseif ii >= 2
tp@0 1358 ivrange = uint32([startindicessinglediff(ii+1):endindicessinglediff(ii+1)]);
tp@0 1359 masterivlistselect = IVNDIFFMATRIX(ivrange,1);
tp@0 1360 masterivlistselect = masterivlistselect(find(singlediffcol(ivrange)==1));
tp@0 1361
tp@0 1362 ivlist = 0;
tp@0 1363 for jj = 3:ii+1
tp@0 1364 ivlist = ivlist + double(POTENTIALISES(masterivlistselect,jj))*ndecimaldivider^(ii+1-jj);
tp@0 1365 end
tp@0 1366
tp@0 1367 IRoriginsfrom(masterivlistselect) = uint32(full(PointertoIRcombs(ivlist)));
tp@0 1368
tp@0 1369 ivlist = ivlist + double(POTENTIALISES(masterivlistselect,2))*ndecimaldivider^(ii-1);
tp@0 1370 A = uint32(ivlist);
tp@0 1371 [B,IA,JA] = unique(A);
tp@0 1372 PointertoIRcombs(B) = masterivlistselect(IA);
tp@0 1373
tp@0 1374 % Now we check if there any active wall numbers that don't
tp@0 1375 % occur as spec1 in an edge-spec1-spec2-spec3 combination.
tp@0 1376
tp@0 1377 wallist = 0;
tp@0 1378 for jj = 1:ii
tp@0 1379 wallist = wallist + double(POTENTIALISES(IVNSPECMATRIX(1:lengthNspecmatrix(ii),ii),jj))*ndecimaldivider^(ii-jj);
tp@0 1380 end
tp@0 1381
tp@0 1382 ivnotincludedyet = find(~ismember(wallist,B));
tp@0 1383
tp@0 1384 if ~isempty(ivnotincludedyet)
tp@0 1385 PointertoIRcombs(wallist(ivnotincludedyet)) = IVNSPECMATRIX(ivnotincludedyet,ii);
tp@0 1386 end
tp@0 1387
tp@0 1388 end
tp@0 1389
tp@0 1390 end
tp@0 1391
tp@0 1392 else
tp@0 1393 ndecimaldivider = 0;
tp@0 1394 PointertoIRcombs = [];
tp@0 1395 IRoriginsfrom = [];
tp@0 1396 end
tp@0 1397
tp@0 1398 %-------------------------------------------------------
tp@0 1399 % Save the output data
tp@0 1400
tp@0 1401 maxval = max(IRoriginsfrom);
tp@0 1402 if maxval < 256
tp@0 1403 IRoriginsfrom = uint8(IRoriginsfrom);
tp@0 1404 elseif maxval < 65536
tp@0 1405 IRoriginsfrom = uint16(IRoriginsfrom);
tp@0 1406 end
tp@0 1407 maxval = max(ORIGINSFROM);
tp@0 1408 if maxval < 256
tp@0 1409 ORIGINSFROM = uint8(ORIGINSFROM);
tp@0 1410 elseif maxval < 65536
tp@0 1411 ORIGINSFROM = uint16(ORIGINSFROM);
tp@0 1412 end
tp@0 1413
tp@0 1414 if SUPPRESSFILES == 0
tp@0 1415 % The complete variable set could be saved if one wanted to.
tp@0 1416 Varlist = [' POTENTIALISES ORIGINSFROM ISCOORDS ISESVISIBILITY IVNSPECMATRIX lengthNspecmatrix IVNDIFFMATRIX lengthNdiffmatrix '];
tp@0 1417 Varlist = [Varlist,' singlediffcol REFLORDER startindicessinglediff endindicessinglediff ndecimaldivider PointertoIRcombs IRoriginsfrom'];
tp@0 1418
tp@0 1419 % Varlist = [' lengthNspecmatrix lengthNdiffmatrix '];
tp@0 1420 % Varlist = [Varlist,' singlediffcol startindicessinglediff endindicessinglediff ndecimaldivider PointertoIRcombs IRoriginsfrom'];
tp@0 1421
tp@0 1422 eval(['save ',ISEStreefile,Varlist])
tp@0 1423 else
tp@0 1424 ISEStreefile = struct('POTENTIALISES',POTENTIALISES,'ORIGINSFROM',ORIGINSFROM,'ISCOORDS',ISCOORDS,'ISESVISIBILITY',ISESVISIBILITY,'IVNSPECMATRIX',IVNSPECMATRIX,...
tp@0 1425 'lengthNspecmatrix',lengthNspecmatrix,'IVNDIFFMATRIX',IVNDIFFMATRIX,'lengthNdiffmatrix',lengthNdiffmatrix,'singlediffcol',singlediffcol,...
tp@0 1426 'REFLORDER',REFLORDER,'startindicessinglediff',startindicessinglediff,'endindicessinglediff',endindicessinglediff,'ndecimaldivider',ndecimaldivider,...
tp@0 1427 'PointertoIRcombs',PointertoIRcombs,'IRoriginsfrom',IRoriginsfrom);
tp@0 1428 end
tp@0 1429
tp@0 1430
tp@0 1431
tp@0 1432