EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_box.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ box geometry headers
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Iwan Kawrakow, 2005
25 #
26 # Contributors: Reid Townson
27 #
28 ###############################################################################
29 */
30 
31 
37 #ifndef EGS_BOX_
38 #define EGS_BOX_
39 
40 #include "egs_base_geometry.h"
41 #include "egs_transformations.h"
42 
43 
44 #ifdef WIN32
45 
46  #ifdef BUILD_BOX_DLL
47  #define EGS_BOX_EXPORT __declspec(dllexport)
48  #else
49  #define EGS_BOX_EXPORT __declspec(dllimport)
50  #endif
51  #define EGS_BOX_LOCAL
52 
53 #else
54 
55  #ifdef HAVE_VISIBILITY
56  #define EGS_BOX_EXPORT __attribute__ ((visibility ("default")))
57  #define EGS_BOX_LOCAL __attribute__ ((visibility ("hidden")))
58  #else
59  #define EGS_BOX_EXPORT
60  #define EGS_BOX_LOCAL
61  #endif
62 
63 #endif
64 
104 class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {
105 
106  EGS_Float ax, ay, az;
107  EGS_AffineTransform *T;
108  static string type;
109 
110 public:
111 
112  EGS_Box(EGS_Float a, const EGS_AffineTransform *t = 0,
113  const string &Name = "") : EGS_BaseGeometry(Name),
114  ax(a), ay(a), az(a), T(0) {
115  if (t) {
116  T = new EGS_AffineTransform(*t);
117  }
118  nreg = 1;
119  };
120 
121  EGS_Box(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
122  const EGS_AffineTransform *t = 0,
123  const string &Name = "") : EGS_BaseGeometry(Name),
124  ax(Ax), ay(Ay), az(Az), T(0) {
125  if (t) {
126  T = new EGS_AffineTransform(*t);
127  }
128  nreg = 1;
129  };
130 
131  ~EGS_Box() {
132  if (T) {
133  delete T;
134  }
135  };
136 
137  bool isInside(const EGS_Vector &x) {
138  EGS_Vector xp = T ? x*(*T) : x;
139  if (2*xp.x + ax < 0 || 2*xp.x - ax > 0) {
140  return false;
141  }
142  if (2*xp.y + ay < 0 || 2*xp.y - ay > 0) {
143  return false;
144  }
145  if (2*xp.z + az < 0 || 2*xp.z - az > 0) {
146  return false;
147  }
148  return true;
149  };
150 
151  int isWhere(const EGS_Vector &x) {
152  return isInside(x) ? 0 : -1;
153  };
154 
155  int inside(const EGS_Vector &x) {
156  return isWhere(x);
157  };
158 
159  EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
160  const EGS_Vector &u) {
161  EGS_Float t = veryFar;
162  howfar(ireg,x,u,t);
163  return t;
164  };
165 
166  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
167  EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
168  EGS_Vector xp(x), up(u);
169  if (T) {
170  xp = x*(*T);
171  up = u*T->getRotation();
172  }
173  if (ireg == 0) {
174  EGS_Float t1 = veryFar;
175  int inew = 0;
176  EGS_Vector n;
177  if (up.x > 0) {
178  t1 = (ax - 2*xp.x)/(2*up.x);
179  n.x = -1;
180  }
181  else if (up.x < 0) {
182  t1 = -(ax + 2*xp.x)/(2*up.x);
183  n.x = 1;
184  }
185  if (t1 < t) {
186  t = t1;
187  inew = -1;
188  }
189  t1 = veryFar;
190  if (up.y > 0) {
191  t1 = (ay - 2*xp.y)/(2*up.y);
192  }
193  else if (up.y < 0) {
194  t1 = -(ay + 2*xp.y)/(2*up.y);
195  }
196  if (t1 < t) {
197  n.x = 0;
198  n.y = up.y > 0 ? -1 : 1;
199  t = t1;
200  inew = -1;
201  }
202  t1 = veryFar;
203  if (up.z > 0) {
204  t1 = (az - 2*xp.z)/(2*up.z);
205  }
206  else if (up.z < 0) {
207  t1 = -(az + 2*xp.z)/(2*up.z);
208  }
209  if (t1 < t) {
210  t = t1;
211  inew = -1;
212  n.x = n.y = 0;
213  n.z = up.z > 0 ? -1 : 1;
214  }
215  if (inew < 0) {
216  if (newmed) {
217  *newmed = -1;
218  }
219  if (normal) {
220  if (T) {
221  *normal = T->getRotation()*n;
222  }
223  else {
224  *normal = n;
225  }
226  }
227  }
228  return inew;
229  }
230  EGS_Float t1 = veryFar;
231  if (2*xp.x + ax < 0 && up.x > 0) {
232  t1 = -(2*xp.x + ax)/(2*up.x);
233  }
234  else if (2*xp.x - ax > 0 && up.x < 0) {
235  t1 = -(2*xp.x - ax)/(2*up.x);
236  }
237  if (t1 < t) {
238  EGS_Float y1 = xp.y + up.y*t1, z1 = xp.z + up.z*t1;
239  if (2*y1 + ay > 0 && 2*y1 - ay < 0 &&
240  2*z1 + az > 0 && 2*z1 - az < 0) {
241  t = t1;
242  if (newmed) {
243  *newmed = med;
244  }
245  if (normal) {
246  if (T) {
247  const EGS_RotationMatrix &R = T->getRotation();
248  if (up.x > 0) {
249  *normal = R*EGS_Vector(-1,0,0);
250  }
251  else {
252  *normal = R*EGS_Vector(1,0,0);
253  }
254  }
255  else {
256  if (up.x > 0) {
257  *normal = EGS_Vector(-1,0,0);
258  }
259  else {
260  *normal = EGS_Vector(1,0,0);
261  }
262  }
263  }
264  return 0;
265  }
266  }
267  t1 = veryFar;
268  if (2*xp.y + ay < 0 && up.y > 0) {
269  t1 = -(2*xp.y + ay)/(2*up.y);
270  }
271  else if (2*xp.y - ay > 0 && up.y < 0) {
272  t1 = -(2*xp.y - ay)/(2*up.y);
273  }
274  if (t1 < t) {
275  EGS_Float x1 = xp.x + up.x*t1, z1 = xp.z + up.z*t1;
276  if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
277  2*z1 + az > 0 && 2*z1 - az < 0) {
278  t = t1;
279  if (newmed) {
280  *newmed = med;
281  }
282  if (normal) {
283  if (T) {
284  const EGS_RotationMatrix &R = T->getRotation();
285  if (up.y > 0) {
286  *normal = R*EGS_Vector(0,-1,0);
287  }
288  else {
289  *normal = R*EGS_Vector(0,1,0);
290  }
291  }
292  else {
293  if (up.y > 0) {
294  *normal = EGS_Vector(0,-1,0);
295  }
296  else {
297  *normal = EGS_Vector(0,1,0);
298  }
299  }
300  }
301  return 0;
302  }
303  }
304  t1 = veryFar;
305  if (2*xp.z + az < 0 && up.z > 0) {
306  t1 = -(2*xp.z + az)/(2*up.z);
307  }
308  else if (2*xp.z - az > 0 && up.z < 0) {
309  t1 = -(2*xp.z - az)/(2*up.z);
310  }
311  if (t1 < t) {
312  EGS_Float x1 = xp.x + up.x*t1, y1 = xp.y + up.y*t1;
313  if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
314  2*y1 + ay > 0 && 2*y1 - ay < 0) {
315  t = t1;
316  if (newmed) {
317  *newmed = med;
318  }
319  if (normal) {
320  if (T) {
321  const EGS_RotationMatrix &R = T->getRotation();
322  if (up.z > 0) {
323  *normal = R*EGS_Vector(0,0,-1);
324  }
325  else {
326  *normal = R*EGS_Vector(0,0,1);
327  }
328  }
329  else {
330  if (up.z > 0) {
331  *normal = EGS_Vector(0,0,-1);
332  }
333  else {
334  *normal = EGS_Vector(0,0,1);
335  }
336  }
337  }
338  return 0;
339  }
340  }
341  return -1;
342  };
343 
344  EGS_Float hownear(int ireg, const EGS_Vector &x) {
345  EGS_Vector xp = T ? x*(*T) : x;
346  if (ireg >= 0) {
347  EGS_Float tmin;
348  EGS_Float t1 = xp.x + 0.5*ax, t2 = 0.5*ax - xp.x;
349  if (t1 < t2) {
350  tmin = t1;
351  }
352  else {
353  tmin = t2;
354  }
355  t1 = xp.y + 0.5*ay;
356  if (t1 < tmin) {
357  tmin = t1;
358  }
359  t1 = 0.5*ay - xp.y;
360  if (t1 < tmin) {
361  tmin = t1;
362  }
363  t1 = xp.z + 0.5*az;
364  if (t1 < tmin) {
365  tmin = t1;
366  }
367  t1 = 0.5*az - xp.z;
368  if (t1 < tmin) {
369  tmin = t1;
370  }
371  return tmin;
372  }
373  EGS_Float s1=0, s2=0;
374  int nout = 0;
375  if (2*xp.x + ax < 0) {
376  EGS_Float t = -0.5*ax - xp.x;
377  s1 += t;
378  s2 += t*t;
379  nout++;
380  }
381  else if (2*xp.x - ax > 0) {
382  EGS_Float t = xp.x - 0.5*ax;
383  s1 += t;
384  s2 += t*t;
385  nout++;
386  }
387  if (2*xp.y + ay < 0) {
388  EGS_Float t = -0.5*ay - xp.y;
389  s1 += t;
390  s2 += t*t;
391  nout++;
392  }
393  else if (2*xp.y - ay > 0) {
394  EGS_Float t = xp.y - 0.5*ay;
395  s1 += t;
396  s2 += t*t;
397  nout++;
398  }
399  if (2*xp.z + az < 0) {
400  EGS_Float t = -0.5*az - xp.z;
401  s1 += t;
402  s2 += t*t;
403  nout++;
404  }
405  else if (2*xp.z - az > 0) {
406  EGS_Float t = xp.z - 0.5*az;
407  s1 += t;
408  s2 += t*t;
409  nout++;
410  }
411  if (nout < 2) {
412  return s1;
413  }
414  return sqrt(s2);
415  };
416 
417  const string &getType() const {
418  return type;
419  };
420 
421  void printInfo() const;
422 
423 };
424 
425 #endif
EGS_BaseGeometry::getType
virtual const string & getType() const =0
Get the geometry type.
EGS_BaseGeometry::isInside
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
EGS_BaseGeometry::printInfo
virtual void printInfo() const
Print information about this geometry.
Definition: egs_base_geometry.cpp:726
EGS_BaseGeometry::howfarToOutside
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
Definition: egs_base_geometry.cpp:381
EGS_AffineTransform
A class providing affine transformations.
Definition: egs_transformations.h:417
EGS_RotationMatrix
A class for vector rotations.
Definition: egs_transformations.h:64
EGS_BaseGeometry::EGS_BaseGeometry
EGS_BaseGeometry(const string &Name)
Construct a geometry named Name.
Definition: egs_base_geometry.cpp:429
EGS_Vector::y
EGS_Float y
y-component
Definition: egs_vector.h:61
EGS_Vector
A class representing 3D vectors.
Definition: egs_vector.h:56
EGS_BaseGeometry::isWhere
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
EGS_BaseGeometry::howfar
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
veryFar
const EGS_Float veryFar
A very large float.
Definition: egs_functions.h:104
EGS_BaseGeometry::med
int med
Medium index.
Definition: egs_base_geometry.h:764
EGS_BaseGeometry
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
Definition: egs_base_geometry.h:93
EGS_Vector::z
EGS_Float z
z-component
Definition: egs_vector.h:62
EGS_BaseGeometry::inside
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
egs_transformations.h
EGS_AffineTransform and EGS_RotationMatrix class header file.
EGS_BaseGeometry::hownear
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
EGS_Vector::x
EGS_Float x
x-component
Definition: egs_vector.h:60
EGS_Box
A box geometry.
Definition: egs_box.h:104
egs_base_geometry.h
EGS_BaseGeometry class header file.
Generated by   doxygen 1.8.5