EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_spherical_shell.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ spherical shell shape
5 # Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
6 # Marc J. P. Chamberland, D. W. O. Rogers
7 #
8 # This file is part of EGSnrc.
9 #
10 # EGSnrc is free software: you can redistribute it and/or modify it under
11 # the terms of the GNU Affero General Public License as published by the
12 # Free Software Foundation, either version 3 of the License, or (at your
13 # option) any later version.
14 #
15 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
16 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
18 # more details.
19 #
20 # You should have received a copy of the GNU Affero General Public License
21 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
22 #
23 ###############################################################################
24 #
25 # Author: Randle Taylor, 2016
26 #
27 # Contributors: Marc Chamberland
28 # Rowan Thomson
29 # Dave Rogers
30 #
31 ###############################################################################
32 #
33 # egs_spherical_shell was developed for the Carleton Laboratory for
34 # Radiotherapy Physics.
35 #
36 ###############################################################################
37 */
38 
39 
45 #include "egs_spherical_shell.h"
46 #include "egs_input.h"
47 #include "egs_functions.h"
48 
49 EGS_SphericalShellShape::EGS_SphericalShellShape(EGS_Float ri, EGS_Float ro, int hemisph, EGS_Float halfangle, const EGS_Vector &Xo,
50  const string &Name,EGS_ObjectFactory *f) :
51  EGS_BaseShape(Name, f), r_inner(ri), r_outer(ro), hemisphere(hemisph), half_angle(halfangle), xo(Xo) {
52  otype = "sphericalShell";
53  if (half_angle < 0) {
54  sgn = -1;
55  half_angle = fabs(half_angle);
56  }
57  else {
58  sgn = 1;
59  }
60 };
61 
62 
65 
66  EGS_Float rnd1 = rndm->getUniform();
67  EGS_Float rnd2 = rndm->getUniform();
68 
69  EGS_Float rad = r_inner + (r_outer - r_inner)*pow(rnd1, 1/3.);
70 
71  EGS_Float cos_max, cos_min;
72  EGS_Float cos_th;
73 
74  if (fabs(half_angle) < 1E-4) {
75  cos_th = 2*rnd2 -1;
76  }
77  else {
78  /* conical section */
79  cos_min = cos(half_angle);
80  cos_th = cos_min + rnd2*(1 - cos_min);
81  rad *= sgn;
82  }
83 
84  EGS_Float sin_th = sqrt(1-cos_th*cos_th);
85 
86  EGS_Float cos_phi, sin_phi;
87  rndm->getAzimuth(cos_phi, sin_phi);
88 
89  EGS_Float x = rad*sin_th*cos_phi;
90  EGS_Float y = rad*sin_th*sin_phi;
91  EGS_Float z = rad*cos_th;
92 
93  if (hemisphere != 0) {
94  z = hemisphere*fabs(z);
95  }
96 
97  return xo + EGS_Vector(x, y, z);
98 
99 };
100 
101 
103  EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
104  EGS_Vector xo = T ? Xo*(*T) : Xo;
105  EGS_Float cost = 2*rndm->getUniform()-1;
106  EGS_Float sint = 1-cost*cost;
107  EGS_Vector x;
108  if (sint > 1e-10) {
109  EGS_Float cphi, sphi;
110  rndm->getAzimuth(cphi,sphi);
111  sint = r_outer*sqrt(sint);
112  x.x = sint*cphi;
113  x.y = sint*sphi;
114  x.z = r_outer*cost;
115  }
116  else {
117  x.z = r_outer*cost;
118  }
119  u = (x + this->xo) - xo;
120  EGS_Float di = 1/u.length();
121  u *= di;
122  wt = u*x*4*M_PI*r_outer*di*di;
123 };
124 
125 EGS_Float EGS_SphericalShellShape::area() const {
126  return 4*M_PI*r_outer*r_outer;
127 };
128 
129 extern "C" {
130 
131  EGS_SPHERICAL_SHELL_EXPORT EGS_BaseShape *createShape(EGS_Input *input, EGS_ObjectFactory *f) {
132 
133  if (!input) {
134  egsWarning("createShape(sphericalShell): null input?\n");
135  return 0;
136  }
137 
138  EGS_Float ri, ro;
139  int err = input->getInput("inner radius", ri);
140  if (err) {
141  egsWarning("createShape(sphericalShell): no 'inner radius' input\n");
142  return 0;
143  }
144  else if (ri < 0) {
145  egsWarning("createShape(sphericalShell): 'inner radius' must be >= 0\n");
146  return 0;
147  }
148 
149  err = input->getInput("outer radius", ro);
150  if (err) {
151  egsWarning("createShape(sphericalShell): no 'outer radius' input\n");
152  return 0;
153  }
154  else if (ro < 0) {
155  egsWarning("createShape(sphericalShell): 'outer radius' must be >= 0\n");
156  return 0;
157  }
158  else if (ri > ro) {
159  egsWarning("createShape(sphericalShell): 'inner radius' must be less than 'outer radius'\n");
160  return 0;
161  }
162 
163  int hemisphere;
164  err = input->getInput("hemisphere", hemisphere);
165  if (err) {
166  hemisphere = 0;
167  }
168  else if ((hemisphere != 0) && (hemisphere != -1) && (hemisphere != 1)) {
169  egsWarning("createShape(sphericalShell): 'hemisphere' must be 1, -1, 0\n");
170  return 0;
171  }
172 
173  EGS_Float half_angle=0, half_angle_deg= 0;
174  err = input->getInput("half angle", half_angle_deg);
175  if (err) {
176  half_angle = 0;
177  }
178  else {
179  half_angle = M_PI/180.*half_angle_deg;
180  }
181 
182  if (half_angle && hemisphere) {
183  egsWarning("createShape(sphericalShell): Hemisphere and half angle specified! Remove one and try again.");
184  return 0;
185  }
186 
187 
188  EGS_SphericalShellShape *result;
189 
190  vector<EGS_Float> xo;
191  err = input->getInput("midpoint",xo);
192  if (err || xo.size() != 3) {
193  xo.clear();
194  xo.push_back(0);
195  xo.push_back(0);
196  xo.push_back(0);
197  }
198 
199  result = new EGS_SphericalShellShape(ri, ro, hemisphere, half_angle, EGS_Vector(xo[0],xo[1],xo[2]));
200  result->setName(input);
201  return result;
202  }
203 
204 }
EGS_Float sgn
The sphere radius.
EGS_Input class header file.
EGS_Float y
y-component
Definition: egs_vector.h:61
A class representing 3D vectors.
Definition: egs_vector.h:56
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
Global egspp functions header file.
Base shape class. All shapes in the EGSnrc C++ class library are derived from EGS_BaseShape.
Definition: egs_shapes.h:112
EGS_SphericalShellShape(EGS_Float ri, EGS_Float ro, int hemisph=0, EGS_Float halfangle=0, 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.
void getPointSourceDirection(const EGS_Vector &Xo, EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt)
Sets the direction u by picking a random point uniformely on the sphere surface.
Base random number generator class. All random number generators should be derived from this class...
Definition: egs_rndm.h:67
EGS_Float z
z-component
Definition: egs_vector.h:62
a spherical shell shape
EGS_Float area() const
Returns the sphere surface area.
A spherical shell shape.
An object factory.
EGS_AffineTransform * T
The affine transformation attached to the shape.
Definition: egs_shapes.h:234
void setName(EGS_Input *inp)
Set the name of the object from the information provided by inp.
EGS_Float getUniform()
Returns a random number uniformly distributed between zero (inclusive) and 1 (exclusive).
Definition: egs_rndm.h:103
EGS_Float x
x-component
Definition: egs_vector.h:60
EGS_Vector xo
The sphere midpoint.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
string otype
The object type.
EGS_Float half_angle
Half angle of conical section.
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
EGS_Vector getPoint(EGS_RandomGenerator *rndm)
Returns a random point within the spherical shell.