EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_iplanes.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ iplanes 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: Frederic Tessier
27 # Reid Townson
28 #
29 ###############################################################################
30 */
31 
32 
38 #ifndef EGS_IPLANES_
39 #define EGS_IPLANES_
40 
41 #include "egs_base_geometry.h"
42 #include "egs_transformations.h"
43 
44 #ifdef WIN32
45 
46  #ifdef BUILD_IPLANES_DLL
47  #define EGS_IPLANES_EXPORT __declspec(dllexport)
48  #else
49  #define EGS_IPLANES_EXPORT __declspec(dllimport)
50  #endif
51  #define EGS_IPLANES_LOCAL
52 
53 #else
54 
55  #ifdef HAVE_VISIBILITY
56  #define EGS_IPLANES_EXPORT __attribute__ ((visibility ("default")))
57  #define EGS_IPLANES_LOCAL __attribute__ ((visibility ("hidden")))
58  #else
59  #define EGS_IPLANES_EXPORT
60  #define EGS_IPLANES_LOCAL
61  #endif
62 
63 #endif
64 
150 class EGS_IPLANES_EXPORT EGS_IPlanes : public EGS_BaseGeometry {
151 
152 public:
153 
161  EGS_IPlanes(const EGS_Vector &Xo, const EGS_Vector &A, int np,
162  const EGS_Float *angles, const string &Name = "",
163  bool degree=true);
164 
172  EGS_IPlanes(const EGS_Vector &Xo, const EGS_Vector &A, int np,
173  const EGS_Vector *aj, const EGS_Float *dj,
174  const string &Name = "");
175 
183  EGS_IPlanes(const EGS_Vector &Xo, const EGS_Vector &A, int np,
184  EGS_Float first=0, const string &Name = "");
185 
186  ~EGS_IPlanes();
187 
188  bool isInside(const EGS_Vector &x) {
189  return true;
190  };
191  int isWhere(const EGS_Vector &x);
192 
193  int inside(const EGS_Vector &x);
194 
195  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
196  EGS_Float &t, int *newmed = 0, EGS_Vector *normal = 0);
197 
198  EGS_Float hownear(int ireg, const EGS_Vector &x);
199 
200  const string &getType() const {
201  return type;
202  };
203 
204  void printInfo() const;
205 
206  EGS_Vector getAxisXo() const {
207  return xo;
208  };
209  EGS_Vector getAxisDirection() const {
210  return axis;
211  };
212  EGS_Vector getPlaneNormal(int j) const {
213  return a[j];
214  };
215  EGS_Float getPlanePosition(int j) const {
216  return d[j];
217  };
218 
219 private:
220 
221  EGS_Vector xo;
222  EGS_Vector axis;
223  EGS_Vector *a;
224  EGS_Float *d;
225  static string type;
226 };
227 
341 class EGS_IPLANES_EXPORT EGS_RadialRepeater : public EGS_BaseGeometry {
342 
343 public:
344 
345  EGS_RadialRepeater(const EGS_Vector &Xo, const EGS_Vector &A, int np,
346  EGS_BaseGeometry *G, EGS_Float first=0, const string &Name = "");
347 
349 
350  bool isInside(const EGS_Vector &x) {
351  return true;
352  };
353 
354  int isWhere(const EGS_Vector &x) {
355  int ir = iplanes->isWhere(x);
356  //EGS_Vector xp = x*R[ir];
357  EGS_Vector xp = (x-xo)*R[ir];
358  int il = g->isWhere(xp);
359  return il >= 0 ? ir*ng + il : nreg-1;
360  };
361 
362  int medium(int ireg) const {
363  if (ireg < 0 || ireg >= nreg) {
364  return -1;
365  }
366  if (ireg == nreg-1) {
367  return med;
368  }
369  int ir = ireg/ng;
370  int il = ireg - ir*ng;
371  return g->medium(il);
372  };
373 
374  int inside(const EGS_Vector &x) {
375  return isWhere(x);
376  };
377 
378  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
379  EGS_Float &t, int *newmed = 0, EGS_Vector *normal = 0) {
380  if (ireg < 0) {
381  egsFatal("\nEGS_RadialRepeater::howfar: position can not"
382  " be outside\n");
383  return ireg;
384  }
385  //egsWarning("\nhowfar(ir=%d x=(%g,%g,%g) u=(%g,%g,%g) t=%g)\n",
386  // ireg,x.x,x.y,x.z,u.x,u.y,u.z,t);
387  if (ireg < nreg-1) { // in the repeated geometry
388  int ir = ireg/ng;
389  int il = ireg - ir*ng;
390  //EGS_Vector xp = x*R[ir], up = u*R[ir];
391  EGS_Vector xp = (x-xo)*R[ir], up = u*R[ir];
392  //egsWarning("In repetition %d: xp=(%g,%g,%g) up=(%g,%g,%g)\n",
393  // ir,xp.x,xp.y,xp.z,up.x,up.y,up.z);
394  int inew = g->howfar(il,xp,up,t,newmed,normal);
395  //egsWarning("il=%d inew=%d t=%d\n",il,inew,t);
396  if (inew < 0) {
397  if (normal) {
398  *normal = R[ir]*(*normal);
399  }
400  inew = nreg-1;
401  if (newmed) {
402  *newmed = med;
403  }
404  }
405  else {
406  inew += ir*ng;
407  }
408  return inew;
409  }
410  // outside of the replicas
411  int ir = iplanes->isWhere(x);
412  //egsWarning("outside repetions, ir=%d\n",ir);
413  EGS_Float t_left = t;
414  EGS_Vector xtmp(x);
415  EGS_Float ttot = 0;
416  //EGS_Vector tmp_n;
417  //EGS_Vector *norm = normal ? &tmp_n : 0;
418  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
419  if (loopCount == loopMax) {
420  egsFatal("EGS_RadialRepeater::howfar: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
421  return -1;
422  }
423  EGS_Float this_t = t_left;
424  //EGS_Vector xp = xtmp*R[ir], up = u*R[ir];
425  EGS_Vector xp = (xtmp-xo)*R[ir], up = u*R[ir];
426  //egsWarning("xtmp=(%g,%g,%g)\n",xtmp.x,xtmp.y,xtmp.z);
427  //egsWarning("xp=(%g,%g,%g) up=(%g,%g,%g)\n",xp.x,xp.y,xp.z,up.x,up.y,up.z);
428  int inew = g->howfar(-1,xp,up,this_t,newmed,normal);
429  //egsWarning("inew=%d t=%g\n",inew,this_t);
430  if (inew >= 0) {
431  t = ttot + this_t;
432  if (normal) {
433  *normal = R[ir]*(*normal);
434  }
435  return ir*ng + inew;
436  }
437  int next_ir = iplanes->howfar(ir,xtmp,u,this_t,0,0);
438  //egsWarning("next sector: %d t=%g\n",next_ir,this_t);
439  if (next_ir == ir) {
440  return ireg;
441  }
442  ttot += this_t;
443  t_left -= this_t;
444  xtmp += u*this_t;
445  ir = next_ir;
446  }
447 
448  return ireg;
449  };
450 
451  EGS_Float hownear(int ireg, const EGS_Vector &x) {
452  if (ireg < 0) {
453  egsFatal("EGS_RadialRepeater::hownear: position can not"
454  " be outside\n");
455  return 0;
456  }
457  if (ireg < nreg-1) { // in the repeated geometry
458  int ir = ireg/ng;
459  int il = ireg - ir*ng;
460  EGS_Vector xp = x*R[ir];
461  return g->hownear(il,xp);
462  }
463  // outside of the replicas
464  int ir = iplanes->isWhere(x);
465  EGS_Vector xp = x*R[ir];
466  return g->hownear(-1,xp);
467  };
468 
469  int getMaxStep() const {
470  return nrep*(g->getMaxStep() + 1);
471  };
472 
473  const string &getType() const {
474  return type;
475  };
476 
477  void printInfo() const;
478 
479  void setRLabels(EGS_Input *input);
480  virtual void getLabelRegions(const string &str, vector<int> &regs);
481 
482 protected:
483 
484  EGS_IPlanes *iplanes;
485  EGS_BaseGeometry *g;
487 
488  EGS_Vector xo;
489 
490  int nrep, ng;
491  EGS_Float phi_o;
492 
493  static string type;
494 
495 
496 };
497 
498 
499 #endif
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:95
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
A class for vector rotations.
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
Definition: egs_vector.h:56
A set of planes intersecting in the same axis.
Definition: egs_iplanes.h:150
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
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.
int med
Medium index.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
EGS_AffineTransform and EGS_RotationMatrix class header file.
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
A radial geometry replicator.
Definition: egs_iplanes.h:341
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
EGS_BaseGeometry class header file.