|
tp@0
|
1 function [validISlist,validIScoords,allreflpoints,listguide,listofreflorder] = EDB1speculISES(eddatafile,...
|
|
tp@0
|
2 S,R,lengthNspecmatrix,specorder,visplanesfromoner)
|
|
tp@0
|
3 % EDB1speculISES - Finds the valid specular reflections by checking visibility and obstruction.
|
|
tp@0
|
4 % Finds the valid specular reflections, given a list of
|
|
tp@0
|
5 % potential IS combinations, by checking visibility and obstruction.
|
|
tp@0
|
6 %
|
|
tp@0
|
7 % Input parameters:
|
|
tp@0
|
8 % eddatafile See a description in EDB1edgeo, EDB1srgeo and EDB1mainISES
|
|
tp@0
|
9 % S
|
|
tp@0
|
10 % R
|
|
tp@0
|
11 % lengthNspecmatrix
|
|
tp@0
|
12 % specorder
|
|
tp@0
|
13 % visplanesfromoner
|
|
tp@0
|
14 % POTENTIALISES (global)
|
|
tp@0
|
15 % ISCOORDS (global)
|
|
tp@0
|
16 % ORIGINSFROM (global)
|
|
tp@0
|
17 % IVNSPECMATRIX (global)
|
|
tp@0
|
18 %
|
|
tp@0
|
19 % GLOBAL parameters:
|
|
tp@0
|
20 % SHOWTEXT JJ JJnumbofchars See EDB1mainISES
|
|
tp@0
|
21 %
|
|
tp@0
|
22 % Output parameters
|
|
tp@0
|
23 % validISlist Matrix [nreflections,specorder] of all reflections, with one row for each reflection.
|
|
tp@0
|
24 % The reflecting plane numbers are given, one for each column.
|
|
tp@0
|
25 % validIScoords Matrix [nreflections,3] of all image source coordinates.
|
|
tp@0
|
26 % allreflpoints Matrix [nreflections,(specorder-1)*3] of all hit points in planes.
|
|
tp@0
|
27 % listguide Matrix [nactiveorders,3] which for each row gives
|
|
tp@0
|
28 % 1. the number of valid reflections for one reflection order
|
|
tp@0
|
29 % 2. the first, and 3. the last
|
|
tp@0
|
30 % row in validISlist etc which contain reflections of that order.
|
|
tp@0
|
31 % Which reflection order each row refers to is given
|
|
tp@0
|
32 % by listofreflorder.
|
|
tp@0
|
33 % listofreflorder List [nactiveorder,3] which for each row gives the specular
|
|
tp@0
|
34 % reflection order that the corresponding row in listguide refers to.
|
|
tp@0
|
35 %
|
|
tp@0
|
36 % Uses functions EDB1chkISvisible EDB1checkobstrpaths
|
|
tp@0
|
37 %
|
|
tp@0
|
38 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
39 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
40 %
|
|
tp@0
|
41 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
42 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
43 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
44 %
|
|
tp@0
|
45 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
46 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
47 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
48 %
|
|
tp@0
|
49 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
50 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
51 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
52 % Peter Svensson (svensson@iet.ntnu.no) 20041109
|
|
tp@0
|
53 %
|
|
tp@0
|
54 % [validISlist,validIScoords,allreflpoints,listguide,listofreflorder] = EDB1speculISES(eddatafile,...
|
|
tp@0
|
55 % S,R,lengthNspecmatrix,specorder,visplanesfromoner);
|
|
tp@0
|
56
|
|
tp@0
|
57 global SHOWTEXT JJ JJnumbofchars
|
|
tp@0
|
58 global POTENTIALISES ISCOORDS ORIGINSFROM IVNSPECMATRIX
|
|
tp@0
|
59
|
|
tp@0
|
60 eval(['load ',eddatafile])
|
|
tp@0
|
61
|
|
tp@0
|
62 [ncorners,slask] = size(corners);
|
|
tp@0
|
63 [nplanes,ncornersperplane] = size(planecorners);
|
|
tp@0
|
64
|
|
tp@0
|
65 [n1,n2] = size(POTENTIALISES);
|
|
tp@0
|
66 if n2 < specorder
|
|
tp@0
|
67 disp([' WARNING: The ISEStree-file does not contain high enough order to support'])
|
|
tp@0
|
68 disp([' spec. refl. order ',int2str(specorder)])
|
|
tp@0
|
69 specorder_maxpossible = n2;
|
|
tp@0
|
70 else
|
|
tp@0
|
71 specorder_maxpossible = specorder;
|
|
tp@0
|
72 end
|
|
tp@0
|
73
|
|
tp@0
|
74 %----------------------------------------------
|
|
tp@0
|
75 % Find all IS that can be propagated to higher order.
|
|
tp@0
|
76 % These must come via planes that S can see, and will be stored in xis1.
|
|
tp@0
|
77
|
|
tp@0
|
78 thinlist = find(planeisthin==1);
|
|
tp@0
|
79
|
|
tp@0
|
80 % ###########################################
|
|
tp@0
|
81 % # #
|
|
tp@0
|
82 % # First order or higher #
|
|
tp@0
|
83 % # #
|
|
tp@0
|
84 % ###########################################
|
|
tp@0
|
85
|
|
tp@0
|
86 validISlist = [];
|
|
tp@0
|
87 validIScoords = [];
|
|
tp@0
|
88 allreflpoints = [];
|
|
tp@0
|
89 listguide = zeros(specorder,3);
|
|
tp@0
|
90 listofreflorder = zeros(specorder,1);
|
|
tp@0
|
91 listguide(1,2) = 1;
|
|
tp@0
|
92
|
|
tp@0
|
93 obstructtestneeded = (sum(canplaneobstruct)~=0);
|
|
tp@0
|
94
|
|
tp@0
|
95 for norder = 1:specorder_maxpossible
|
|
tp@0
|
96
|
|
tp@0
|
97 if SHOWTEXT >= 3
|
|
tp@0
|
98 disp([' Reflection order ',int2str(norder)])
|
|
tp@0
|
99 end
|
|
tp@0
|
100
|
|
tp@0
|
101 % Start with all the theoretically possible IS
|
|
tp@0
|
102 masterivlist = IVNSPECMATRIX(1:lengthNspecmatrix(norder),norder);
|
|
tp@0
|
103 possiblecombs = POTENTIALISES(masterivlist,1:norder);
|
|
tp@0
|
104
|
|
tp@0
|
105 % Pick out only those IS that come via a last refl plane that the
|
|
tp@0
|
106 % receiver is in front of
|
|
tp@0
|
107 okcombs = find(visplanesfromoner(possiblecombs(:,norder))==2);
|
|
tp@0
|
108
|
|
tp@0
|
109 masterivlist = masterivlist(okcombs);
|
|
tp@0
|
110 possiblecombs = POTENTIALISES(masterivlist,1:norder);
|
|
tp@0
|
111 possibleIS = ISCOORDS(masterivlist,:);
|
|
tp@0
|
112 nis = length(masterivlist);
|
|
tp@0
|
113
|
|
tp@0
|
114 if SHOWTEXT >= 3
|
|
tp@0
|
115 disp([' ',int2str(nis),' IS found initially for order ',int2str(norder)])
|
|
tp@0
|
116 if SHOWTEXT >= 5
|
|
tp@0
|
117 possiblecombs
|
|
tp@0
|
118 end
|
|
tp@0
|
119 end
|
|
tp@0
|
120
|
|
tp@0
|
121 % Go through iteratively the reflection points, and check if they are
|
|
tp@0
|
122 % inside the reflection plane. Start with the last reflection plane and
|
|
tp@0
|
123 % work backwards.
|
|
tp@0
|
124 for jj = norder:-1:1
|
|
tp@0
|
125 if nis > 0
|
|
tp@0
|
126
|
|
tp@0
|
127 if jj == norder
|
|
tp@0
|
128 fromcoords = possibleIS;
|
|
tp@0
|
129 tocoords = R;
|
|
tp@0
|
130 else
|
|
tp@0
|
131
|
|
tp@0
|
132 eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])
|
|
tp@0
|
133 ivreflist = ORIGINSFROM(masterivlist);
|
|
tp@0
|
134 for kk = jj:norder-2
|
|
tp@0
|
135 ivreflist = ORIGINSFROM(ivreflist);
|
|
tp@0
|
136 end
|
|
tp@0
|
137 fromcoords = ISCOORDS(ivreflist,:);
|
|
tp@0
|
138 end
|
|
tp@0
|
139
|
|
tp@0
|
140 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(possiblecombs(:,jj),4),planenvecs(possiblecombs(:,jj),:),minvals(possiblecombs(:,jj),:),...
|
|
tp@0
|
141 maxvals(possiblecombs(:,jj),:),planecorners(possiblecombs(:,jj),:),corners,ncornersperplanevec(possiblecombs(:,jj)));
|
|
tp@0
|
142 if ~isempty(edgehits) | ~isempty(cornerhits)
|
|
tp@0
|
143 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
|
|
tp@0
|
144 disp(' handled correctly yet.')
|
|
tp@0
|
145 end
|
|
tp@0
|
146 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
|
|
tp@0
|
147
|
|
tp@0
|
148 masterivlist = masterivlist(hitplanes);
|
|
tp@0
|
149 possiblecombs = POTENTIALISES(masterivlist,1:norder);
|
|
tp@0
|
150 possibleIS = ISCOORDS(masterivlist,:);
|
|
tp@0
|
151 if jj < norder
|
|
tp@0
|
152 for kk = jj+1:norder
|
|
tp@0
|
153 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);'])
|
|
tp@0
|
154 end
|
|
tp@0
|
155 end
|
|
tp@0
|
156 nis = length(masterivlist);
|
|
tp@0
|
157
|
|
tp@0
|
158 if SHOWTEXT >= 3
|
|
tp@0
|
159 disp([' ',int2str(nis),' IS survived the visibility test in refl plane number ',int2str(jj)])
|
|
tp@0
|
160 if SHOWTEXT >= 5
|
|
tp@0
|
161 possiblecombs
|
|
tp@0
|
162 end
|
|
tp@0
|
163 end
|
|
tp@0
|
164 end
|
|
tp@0
|
165
|
|
tp@0
|
166 end
|
|
tp@0
|
167
|
|
tp@0
|
168 if obstructtestneeded & nis > 0
|
|
tp@0
|
169 % Check obstructions for all the paths: S -> plane1 -> plane2 -> ...
|
|
tp@0
|
170 % -> planeN -> R
|
|
tp@0
|
171
|
|
tp@0
|
172 for jj = 1:norder+1
|
|
tp@0
|
173
|
|
tp@0
|
174 if nis > 0
|
|
tp@0
|
175 if jj==1
|
|
tp@0
|
176 fromcoords = S;
|
|
tp@0
|
177 startplanes = [];
|
|
tp@0
|
178 else
|
|
tp@0
|
179 startplanes = possiblecombs(:,jj-1);
|
|
tp@0
|
180 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
|
|
tp@0
|
181 end
|
|
tp@0
|
182 if jj == norder+1
|
|
tp@0
|
183 tocoords = R;
|
|
tp@0
|
184 endplanes = [];
|
|
tp@0
|
185 else
|
|
tp@0
|
186 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
|
|
tp@0
|
187 endplanes = possiblecombs(:,jj);
|
|
tp@0
|
188 end
|
|
tp@0
|
189 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
|
|
tp@0
|
190 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
|
|
tp@0
|
191
|
|
tp@0
|
192 if nobstructions > 0
|
|
tp@0
|
193 masterivlist = masterivlist(nonobstructedpaths);
|
|
tp@0
|
194 possiblecombs = POTENTIALISES(masterivlist,1:norder);
|
|
tp@0
|
195 possibleIS = ISCOORDS(masterivlist,:);
|
|
tp@0
|
196 for kk = 1:norder
|
|
tp@0
|
197 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])
|
|
tp@0
|
198 end
|
|
tp@0
|
199 nis = length(masterivlist);
|
|
tp@0
|
200 end
|
|
tp@0
|
201 if SHOWTEXT >= 3
|
|
tp@0
|
202 disp([' ',int2str(nis),' IS survived the obstruction test for path ',int2str(jj)])
|
|
tp@0
|
203 if SHOWTEXT >= 5
|
|
tp@0
|
204 possiblecombs
|
|
tp@0
|
205 end
|
|
tp@0
|
206 end
|
|
tp@0
|
207 end
|
|
tp@0
|
208
|
|
tp@0
|
209 end
|
|
tp@0
|
210
|
|
tp@0
|
211 end
|
|
tp@0
|
212 if nis > 0
|
|
tp@0
|
213 if ~isempty(validISlist)
|
|
tp@0
|
214 [n1,n2] = size(validISlist);
|
|
tp@0
|
215 [n3,n4] = size(possiblecombs);
|
|
tp@0
|
216 validISlist = [[validISlist zeros(n1,n4-n2)];(possiblecombs)];
|
|
tp@0
|
217 else
|
|
tp@0
|
218 validISlist = possiblecombs;
|
|
tp@0
|
219 end
|
|
tp@0
|
220 validIScoords = [validIScoords;possibleIS];
|
|
tp@0
|
221
|
|
tp@0
|
222 newestreflpoints = [];
|
|
tp@0
|
223 for kk = 1:norder
|
|
tp@0
|
224 eval(['newestreflpoints = [newestreflpoints reflpoints',JJ(kk,1:JJnumbofchars(kk)),'];'])
|
|
tp@0
|
225 end
|
|
tp@0
|
226
|
|
tp@0
|
227 if ~isempty(allreflpoints)
|
|
tp@0
|
228 [n1,n2] = size(allreflpoints);
|
|
tp@0
|
229 [n3,n4] = size(newestreflpoints);
|
|
tp@0
|
230 allreflpoints = [[allreflpoints zeros(n1,n4-n2)];newestreflpoints];
|
|
tp@0
|
231 else
|
|
tp@0
|
232 allreflpoints = newestreflpoints;
|
|
tp@0
|
233 end
|
|
tp@0
|
234 end
|
|
tp@0
|
235
|
|
tp@0
|
236 listguide(norder,1) = nis;
|
|
tp@0
|
237 listguide(norder,3) = listguide(norder,2)+nis-1;
|
|
tp@0
|
238 if norder < specorder
|
|
tp@0
|
239 listguide(norder+1,2) = listguide(norder,3)+1;
|
|
tp@0
|
240 end
|
|
tp@0
|
241 listofreflorder(norder) = norder;
|
|
tp@0
|
242
|
|
tp@0
|
243 end
|
|
tp@0
|
244
|
|
tp@0
|
245 iv = find(listguide(:,3)==0);
|
|
tp@0
|
246 listguide(iv,2) = zeros(size(iv));
|
|
tp@0
|
247
|
|
tp@0
|
248 [n1,n2] = size(validISlist);
|
|
tp@0
|
249 if n2 < specorder
|
|
tp@0
|
250 validISlist = [validISlist zeros(n1,specorder-n2)];
|
|
tp@0
|
251 end
|
|
tp@0
|
252
|
|
tp@0
|
253 [n1,n2] = size(allreflpoints);
|
|
tp@0
|
254 if n2 < specorder*3
|
|
tp@0
|
255 allreflpoints = [allreflpoints zeros(n1,specorder*3-n2)];
|
|
tp@0
|
256 end
|
|
tp@0
|
257
|
|
tp@0
|
258 iv = find(listguide(:,1)==0);
|
|
tp@0
|
259 listguide(iv,:) = [];
|
|
tp@0
|
260 listofreflorder(iv) = [];
|
