EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_shapes.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ shapes 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 # Marc Chamberland
28 # Reid Townson
29 #
30 ###############################################################################
31 */
32 
33 
39 #ifndef EGS_SHAPES_
40 #define EGS_SHAPES_
41 
42 #include "egs_vector.h"
43 #include "egs_transformations.h"
44 #include "egs_rndm.h"
45 #include "egs_object_factory.h"
46 
47 #include <string>
48 using std::string;
49 
50 class EGS_Input;
51 
115 
116 public:
117 
119  EGS_BaseShape(const string &Name="",EGS_ObjectFactory *f=0) :
120  EGS_Object(Name,f), T(0) {
121  otype = "base_shape";
122  };
124  virtual ~EGS_BaseShape() {
125  if (T) {
126  delete T;
127  }
128  };
129 
137  if (T) {
138  return (*T)*getPoint(rndm);
139  }
140  else {
141  return getPoint(rndm);
142  }
143  };
144 
152  egsFatal("You need to implement the getPoint function in your "
153  "derived class\n");
154  return EGS_Vector();
155  };
156 
165  void setTransformation(EGS_Input *inp);
166 
172  if (T) {
173  delete T;
174  }
175  T = new EGS_AffineTransform(*t);
176  };
177 
181  return T;
182  };
183 
191  static EGS_BaseShape *createShape(EGS_Input *inp);
192 
199  static EGS_BaseShape *getShape(const string &Name);
200 
206  virtual bool supportsDirectionMethod() const {
207  return false;
208  };
209 
210  /* getNextShapePosition is the equivalent of getNextParticle but for a shape object. Its goal is to determine the next state of the geometry, either by synchronizing itself to the
211  * source time parameter, or by sampling it's own time parameter and updating itself accordingly if the source has provided no time index.
212  *
213  *This function has a non-empty implementation in 2 cases.
214  *1) it is re implemented in any composite shape, where it will call getNextShapePosition on all of its components
215  *2) it is re implemented in the dynamic shape class. This is where the code will find the current (non static) state of the shape. */
216  virtual void getNextShapePosition(EGS_RandomGenerator *rndm) {};
217 
230  virtual void getPointSourceDirection(const EGS_Vector &xo,
231  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
232  egsFatal("getPointSourceDirection: you have to implement this "
233  "method for the %s shape if you want to use it\n",otype.c_str());
234  };
235 
242  virtual EGS_Float area() const {
243  return 1;
244  };
245 
248  virtual void updatePosition(EGS_Float time) { };
249 
250 protected:
251 
253 
254 };
255 
267 
268 public:
269 
271  EGS_SurfaceShape(const string &Name="",EGS_ObjectFactory *f=0) :
272  EGS_BaseShape(Name,f), A(1) {};
278  bool supportsDirectionMethod() const {
279  return true;
280  };
282  EGS_Float area() const {
283  return A;
284  };
290  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
291  EGS_Vector xo = T ? Xo*(*T) : Xo;
292  EGS_Vector x = getPoint(rndm);
293  u = x - xo;
294  EGS_Float d2i = 1/u.length2(), di = sqrt(d2i);
295  u *= di;
296  wt = A*fabs(u.z)*d2i;
297  if (T) {
298  T->rotate(u);
299  }
300  };
301 
302 protected:
303 
309  EGS_Float A;
310 
311 };
312 
329 
330 public:
331 
334  const string &Name="",EGS_ObjectFactory *f=0) :
335  EGS_BaseShape(Name,f), xo(Xo) {
336  otype = "point";
337  };
338  ~EGS_PointShape() { };
341  return xo;
342  };
347 
348 protected:
349 
351 
352 };
353 
371 
372 protected:
373 
374  EGS_Float ax, ay, az;
375 
376 public:
377 
379  EGS_BoxShape(const string &Name="",EGS_ObjectFactory *f=0) :
380  EGS_BaseShape(Name,f), ax(1), ay(1), az(1) {
381  otype="box";
382  };
384  EGS_BoxShape(EGS_Float A, const EGS_AffineTransform *t = 0,
385  const string &Name="",EGS_ObjectFactory *f=0) :
386  EGS_BaseShape(Name,f), ax(A), ay(A), az(A) {
387  if (t) {
388  T = new EGS_AffineTransform(*t);
389  }
390  otype="box";
391  };
393  EGS_BoxShape(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
394  const EGS_AffineTransform *t = 0,
395  const string &Name="",EGS_ObjectFactory *f=0) :
396  EGS_BaseShape(Name,f), ax(Ax), ay(Ay), az(Az) {
397  if (t) {
398  T = new EGS_AffineTransform(*t);
399  }
400  otype="box";
401  };
404 
407  EGS_Vector v(ax*(rndm->getUniform()-0.5),
408  ay*(rndm->getUniform()-0.5),
409  az*(rndm->getUniform()-0.5));
410  return v;
411  };
412 
417 
421  bool supportsDirectionMethod() const {
422  return true;
423  };
424 
430  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
431  EGS_Vector xo = T ? Xo*(*T) : Xo;
432  EGS_Float eta = rndm->getUniform()*area();
433  if (eta < 2*ax*ay) {
434  u.x = ax*(rndm->getUniform()-0.5);
435  u.y = ay*(rndm->getUniform()-0.5);
436  if (eta < ax*ay) {
437  u.z = az/2;
438  wt = u.z - xo.z;
439  }
440  else {
441  u.z = -az/2;
442  wt = xo.z - u.z;
443  }
444  }
445  else if (eta < 2*(ax*ay + ax*az)) {
446  u.x = ax*(rndm->getUniform()-0.5);
447  u.z = az*(rndm->getUniform()-0.5);
448  if (eta < 2*ax*ay + ax*az) {
449  u.y = ay/2;
450  wt = u.y - xo.y;
451  }
452  else {
453  u.y = -ay/2;
454  wt = xo.y - u.y;
455  }
456  }
457  else {
458  eta -= 2*(ax*ay + ax*az);
459  u.y = ay*(rndm->getUniform()-0.5);
460  u.z = az*(rndm->getUniform()-0.5);
461  if (eta < ay*az) {
462  u.x = ax/2;
463  wt = u.x - xo.x;
464  }
465  else {
466  u.x = -ax/2;
467  wt = xo.x - u.x;
468  }
469  }
470  u -= xo;
471  EGS_Float d2 = u.length2(), d = sqrt(d2);
472  u *= (1/d);
473  wt *= (area()/(d2*d));
474  if (T) {
475  T->rotate(u);
476  }
477  };
479  EGS_Float area() const {
480  return 2*(ax*ay + ax*az + ay*az);
481  };
482 
483 };
484 
503 
504 protected:
505 
506  EGS_Float R;
508 
509 public:
510 
512  EGS_SphereShape(const string &Name="",EGS_ObjectFactory *f=0) :
513  EGS_BaseShape(Name,f), R(1), xo() {
514  otype="sphere";
515  };
517  EGS_SphereShape(EGS_Float r, const EGS_Vector &Xo = EGS_Vector(0,0,0),
518  const string &Name="",EGS_ObjectFactory *f=0) :
519  EGS_BaseShape(Name,f), R(r), xo(Xo) {
520  otype = "sphere";
521  };
524 
527  EGS_Float r = rndm->getUniform(), r1 = rndm->getUniform(),
528  r2 = rndm->getUniform();
529  if (r1 > r) {
530  r = r1;
531  }
532  if (r2 > r) {
533  r = r2;
534  }
535  EGS_Float cost = 2*rndm->getUniform()-1;
536  EGS_Float sint = sqrt(1-cost*cost);
537  r1 = R*r*sint;
538  EGS_Float cphi, sphi;
539  rndm->getAzimuth(cphi,sphi);
540  return xo + EGS_Vector(r1*cphi,r1*sphi,R*r*cost);
541  };
542 
547 
551  bool supportsDirectionMethod() const {
552  return true;
553  };
554 
560  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
561  EGS_Vector xo = T ? Xo*(*T) : Xo;
562  EGS_Float cost = 2*rndm->getUniform()-1;
563  EGS_Float sint = 1-cost*cost;
564  EGS_Vector x;
565  if (sint > epsilon) {
566  EGS_Float cphi, sphi;
567  rndm->getAzimuth(cphi,sphi);
568  sint = R*sqrt(sint);
569  x.x = sint*cphi;
570  x.y = sint*sphi;
571  x.z = R*cost;
572  }
573  else {
574  x.z = R*cost;
575  }
576  u = (x + this->xo) - xo;
577  EGS_Float di = 1/u.length();
578  u *= di;
579  wt = u*x*4*M_PI*R*di*di;
580  };
582  EGS_Float area() const {
583  return 4*M_PI*R*R;
584  };
585 };
586 
611 
612 protected:
613 
614  EGS_Float R;
615  EGS_Float h;
618  EGS_Float phi_min;
619  EGS_Float phi_max;
620  bool has_phi;
621 
623  inline void getPointInCircle(EGS_RandomGenerator *rndm, EGS_Float &x,
624  EGS_Float &y) {
625  if (!has_phi) {
626  do {
627  x = 2*rndm->getUniform()-1;
628  y = 2*rndm->getUniform()-1;
629  }
630  while (x*x + y*y > 1);
631  x *= R;
632  y *= R;
633  }
634  else {
635  EGS_Float r = R*sqrt(rndm->getUniform());
636  EGS_Float eta = rndm->getUniform();
637  EGS_Float phi = phi_min*(1-eta) + phi_max*eta;
638  x = r*cos(phi);
639  y = r*sin(phi);
640  }
641  };
642 
643 
644 public:
645 
649  EGS_CylinderShape(const string &Name="",EGS_ObjectFactory *f=0) :
650  EGS_BaseShape(), R(1), h(1), xo(), a(0,0,1),
651  phi_min(0), phi_max(2*M_PI), has_phi(false) {
652  otype="cylinder";
653  };
654  EGS_CylinderShape(EGS_Float r, EGS_Float H,
655  const EGS_Vector &Xo = EGS_Vector(0,0,0),
656  const EGS_Vector &A = EGS_Vector(0,0,1),
657  const string &Name="",EGS_ObjectFactory *f=0) :
658  EGS_BaseShape(Name,f), R(r), h(H), xo(Xo), a(A),
659  phi_min(0), phi_max(2*M_PI), has_phi(false) {
660  EGS_RotationMatrix rmat(a);
661  if (xo.length2() > epsilon || !rmat.isI()) {
662  T = new EGS_AffineTransform(rmat.inverse(),xo);
663  }
664  otype="cylinder";
665  };
669  EGS_CylinderShape(EGS_Float r, EGS_Float H, const EGS_AffineTransform *t,
670  const string &Name="",EGS_ObjectFactory *f=0) :
671  EGS_BaseShape(Name,f), R(r), h(H),
672  phi_min(0), phi_max(2*M_PI), has_phi(false) {
673  if (t) {
674  T = new EGS_AffineTransform(*t);
675  }
676  otype="cylinder";
677  };
680 
682  void setPhiRange(EGS_Float Phi_min, EGS_Float Phi_max) {
683  if (Phi_min < Phi_max) {
684  phi_min = Phi_min;
685  phi_max = Phi_max;
686  }
687  else {
688  phi_min = Phi_max;
689  phi_max = Phi_min;
690  }
691  if (phi_max - phi_min < 1.99999*M_PI) {
692  has_phi = true;
693  }
694  else {
695  has_phi = false;
696  }
697  };
698 
703  EGS_Float x,y;
704  getPointInCircle(rndm,x,y);
705  EGS_Float z = h*(rndm->getUniform()-0.5);
706  return EGS_Vector(x,y,z);
707  };
708 
713 
715  EGS_Float getRadius() const {
716  return R;
717  };
719  EGS_Float getHeight() const {
720  return h;
721  };
722 
726  bool supportsDirectionMethod() const {
727  return true;
728  };
729 
735  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
736  EGS_Vector xo = T ? Xo*(*T) : Xo;
737  EGS_Float eta = rndm->getUniform()*(R+h);
738  EGS_Vector x; // point on cylinder with respect to midpoint
739  EGS_Vector n; // normal to cylinder at point x
740  if (eta < R) {
741  getPointInCircle(rndm,x.x,x.y);
742  if (2*eta < R) {
743  x.z = h/2; // top face: normal is up
744  n.z = 1;
745  }
746  else {
747  x.z = -h/2; // bottom face: normal is down
748  n.z = -1;
749  }
750  }
751  else {
752  EGS_Float cphi,sphi;
753  rndm->getAzimuth(cphi,sphi);
754  x.x = R*cphi;
755  x.y = R*sphi;
756  x.z = h*(rndm->getUniform()-0.5);
757  n.x = x.x;
758  n.y = x.y; // side face: normal is (x,y)
759  }
760  u = (x+this->xo) - xo; // direction vector from origin to cylinder point
761  EGS_Float d2 = u.length2(), d = sqrt(d2);
762  u *= (1/d); // normalize direction vectors
763  n.normalize(); // normalize normal
764  wt = u*n*area()/d2;
765  if (T) {
766  T->rotate(u);
767  }
768  };
770  EGS_Float area() const {
771  return 2*M_PI*R*(R+h);
772  };
773 };
774 
775 #endif
A class providing affine transformations.
Base shape class. All shapes in the EGSnrc C++ class library are derived from EGS_BaseShape.
Definition: egs_shapes.h:114
void setTransformation(EGS_AffineTransform *t)
Set the transformation attached to this shape.
Definition: egs_shapes.h:171
virtual void getPointSourceDirection(const EGS_Vector &xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Definition: egs_shapes.h:230
virtual EGS_Float area() const
Definition: egs_shapes.h:242
EGS_AffineTransform * T
The affine transformation attached to the shape.
Definition: egs_shapes.h:248
EGS_BaseShape(const string &Name="", EGS_ObjectFactory *f=0)
Construct a shape named Name.
Definition: egs_shapes.h:119
virtual ~EGS_BaseShape()
Destructor. Deletes T if it is not null.
Definition: egs_shapes.h:124
virtual EGS_Vector getPoint(EGS_RandomGenerator *rndm)
Sample and return a random 3D vector.
Definition: egs_shapes.h:151
virtual void updatePosition(EGS_Float time)
Update the position of the shape if it is in motion.
Definition: egs_shapes.h:248
const EGS_AffineTransform * getTransform() const
Get a pointer to the affine transformation attached to this shape.
Definition: egs_shapes.h:180
virtual EGS_Vector getRandomPoint(EGS_RandomGenerator *rndm)
Returns a random 3D vector.
Definition: egs_shapes.h:136
virtual bool supportsDirectionMethod() const
Definition: egs_shapes.h:206
A box shape.
Definition: egs_shapes.h:370
void getPointSourceDirection(const EGS_Vector &Xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Sets the direction u by picking a random point uniformly the on the box surface.
Definition: egs_shapes.h:429
EGS_Float az
The box size.
Definition: egs_shapes.h:374
EGS_Vector getPoint(EGS_RandomGenerator *rndm)
Returns a point uniformly distributed within the box.
Definition: egs_shapes.h:406
EGS_BoxShape(const string &Name="", EGS_ObjectFactory *f=0)
Create a box shape with unit size.
Definition: egs_shapes.h:379
EGS_BoxShape(EGS_Float A, const EGS_AffineTransform *t=0, const string &Name="", EGS_ObjectFactory *f=0)
Create a cube with size A.
Definition: egs_shapes.h:384
EGS_BoxShape(EGS_Float Ax, EGS_Float Ay, EGS_Float Az, const EGS_AffineTransform *t=0, const string &Name="", EGS_ObjectFactory *f=0)
Create a box shape with size Ax,Ay,Az.
Definition: egs_shapes.h:393
bool supportsDirectionMethod() const
Returns true. (It is easy to implement the getPointSourceDirection() method for a box....
Definition: egs_shapes.h:421
EGS_Float area() const
Returns the box surface area.
Definition: egs_shapes.h:479
~EGS_BoxShape()
Destructor. Does nothing.
Definition: egs_shapes.h:403
A cylinder shape.
Definition: egs_shapes.h:610
void getPointSourceDirection(const EGS_Vector &Xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Sets the direction u by picking a random point uniformly on the cylinder surface.
Definition: egs_shapes.h:734
bool supportsDirectionMethod() const
Returns true. (It is easy to implement the getPointSourceDirection() method for a cylinder....
Definition: egs_shapes.h:726
EGS_Float getHeight() const
Definition: egs_shapes.h:719
void getPointInCircle(EGS_RandomGenerator *rndm, EGS_Float &x, EGS_Float &y)
Get a point uniformly distributed within a circle.
Definition: egs_shapes.h:623
void setPhiRange(EGS_Float Phi_min, EGS_Float Phi_max)
Definition: egs_shapes.h:682
EGS_Float h
Cylinder height.
Definition: egs_shapes.h:615
EGS_Float getRadius() const
Definition: egs_shapes.h:715
EGS_CylinderShape(EGS_Float r, EGS_Float H, const EGS_AffineTransform *t, const string &Name="", EGS_ObjectFactory *f=0)
Definition: egs_shapes.h:669
EGS_Vector a
Cylinder axis.
Definition: egs_shapes.h:617
EGS_CylinderShape(const string &Name="", EGS_ObjectFactory *f=0)
Definition: egs_shapes.h:649
bool has_phi
True, if azimuthal range restricted.
Definition: egs_shapes.h:620
EGS_Vector getPoint(EGS_RandomGenerator *rndm)
Samples and returns a point uniformly distributed within the cylinder.
Definition: egs_shapes.h:702
EGS_Float R
Cylinder radius.
Definition: egs_shapes.h:614
EGS_Vector xo
midpoint
Definition: egs_shapes.h:616
EGS_Float area() const
Returns the cylinder surface area.
Definition: egs_shapes.h:770
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
An object factory.
Base egspp object.
virtual EGS_Object * createObject(EGS_Input *inp)
Create an object from the infromation pointed to by inp.
A point shape. This is the simplest shape possible: it simply always returns the same point.
Definition: egs_shapes.h:328
EGS_Vector xo
The point position.
Definition: egs_shapes.h:350
EGS_Vector getPoint(EGS_RandomGenerator *)
Returns a fixed point.
Definition: egs_shapes.h:340
EGS_PointShape(const EGS_Vector &Xo=EGS_Vector(), const string &Name="", EGS_ObjectFactory *f=0)
Construct a point shape located at Xo.
Definition: egs_shapes.h:333
Base random number generator class. All random number generators should be derived from this class.
Definition: egs_rndm.h:67
void getAzimuth(EGS_Float &cphi, EGS_Float &sphi)
Sets cphi and sphi to the cosine and sine of a random angle uniformely distributed between 0 and .
Definition: egs_rndm.h:138
EGS_Float getUniform()
Returns a random number uniformly distributed between zero (inclusive) and 1 (exclusive).
Definition: egs_rndm.h:103
A class for vector rotations.
A sphere shape.
Definition: egs_shapes.h:502
EGS_SphereShape(const string &Name="", EGS_ObjectFactory *f=0)
Construct a sphere of unit radius about the origin.
Definition: egs_shapes.h:512
EGS_SphereShape(EGS_Float r, const EGS_Vector &Xo=EGS_Vector(0, 0, 0), const string &Name="", EGS_ObjectFactory *f=0)
Construct a sphere of radius r with midpoint Xo.
Definition: egs_shapes.h:517
bool supportsDirectionMethod() const
Returns true. (It is easy to implement the getPointSourceDirection() method for a sphere....
Definition: egs_shapes.h:551
void getPointSourceDirection(const EGS_Vector &Xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Sets the direction u by picking a random point uniformly on the sphere surface.
Definition: egs_shapes.h:559
EGS_Float area() const
Returns the sphere surface area.
Definition: egs_shapes.h:582
EGS_Vector getPoint(EGS_RandomGenerator *rndm)
Returns a random point within the sphere.
Definition: egs_shapes.h:526
EGS_Float R
The sphere radius.
Definition: egs_shapes.h:506
EGS_Vector xo
The sphere midpoint.
Definition: egs_shapes.h:507
A surface shape.
Definition: egs_shapes.h:266
EGS_SurfaceShape(const string &Name="", EGS_ObjectFactory *f=0)
Construct a surface shape named Name.
Definition: egs_shapes.h:271
~EGS_SurfaceShape()
Destructor. Does nothing.
Definition: egs_shapes.h:274
void getPointSourceDirection(const EGS_Vector &Xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Get a random direction given a source position Xo.
Definition: egs_shapes.h:289
EGS_Float area() const
Returns the area of this surface shape.
Definition: egs_shapes.h:282
bool supportsDirectionMethod() const
Always returns true. Shapes derived from this class must implement the getPoint() method to return po...
Definition: egs_shapes.h:278
A class representing 3D vectors.
Definition: egs_vector.h:57
EGS_Float y
y-component
Definition: egs_vector.h:62
EGS_Float z
z-component
Definition: egs_vector.h:63
EGS_Float x
x-component
Definition: egs_vector.h:61
#define EGS_EXPORT
Export symbols from the egspp library.
Definition: egs_libconfig.h:90
EGS_Object and EGS_ObjectFactory class header file.
EGS_RandomGenerator class header file.
EGS_AffineTransform and EGS_RotationMatrix class header file.
EGS_Vector methods for the manipulation of 3D vectors in cartesian co-ordinates.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:62