EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_iplanes.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ iplanes geometry
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 #
28 ###############################################################################
29 */
30 
31 
37 #include "egs_iplanes.h"
38 #include "egs_input.h"
39 #include "egs_transformations.h"
40 #include "egs_functions.h"
41 #include "egs_math.h"
42 
43 #include <vector>
44 
45 using namespace std;
46 
47 string EGS_IPlanes::type = "EGS_IPlanes";
48 string EGS_RadialRepeater::type = "EGS_RadialRepeater";
49 
50 EGS_IPlanes::EGS_IPlanes(const EGS_Vector &Xo, const EGS_Vector &A, int np,
51  const EGS_Float *angles, const string &Name, bool degree) :
52  EGS_BaseGeometry(Name), xo(Xo), axis(A) {
53  nreg = 2*np;
54  a = new EGS_Vector[nreg];
55  d = new EGS_Float[nreg];
56  int j;
57  EGS_Float phi180 = M_PI;
58  if (degree) {
59  phi180 = 180;
60  }
61  EGS_RotationMatrix R(axis);
62  for (j=0; j<nreg; j++) {
63  EGS_Float phi;
64  if (j < np) {
65  phi = angles[j];
66  }
67  else {
68  phi = angles[j-np] + phi180;
69  }
70  if (degree) {
71  phi *= (M_PI/180);
72  }
73  EGS_Float cphi = cos(phi), sphi = sin(phi);
74  a[j] = EGS_Vector(-sphi,cphi,0)*R;
75  d[j] = a[j]*xo;
76  }
77 }
78 
80  int np, const EGS_Vector *aj, const EGS_Float *dj,
81  const string &Name) : EGS_BaseGeometry(Name), xo(Xo), axis(A) {
82  nreg = np;
83  a = new EGS_Vector[nreg];
84  d = new EGS_Float[nreg];
85  int j;
86  for (j=0; j<nreg; j++) {
87  a[j] = aj[j];
88  d[j] = dj[j];
89  }
90 }
91 
93  int np, EGS_Float first, const string &Name) :
94  EGS_BaseGeometry(Name), xo(Xo), axis(A) {
95  nreg = np;
96  EGS_Float dphi = 2*M_PI/nreg;
97  EGS_RotationMatrix R(axis);
98  a = new EGS_Vector[nreg];
99  d = new EGS_Float[nreg];
100  for (int j=0; j<nreg; j++) {
101  EGS_Float phi = first + dphi*j;
102  EGS_Float cphi = cos(phi), sphi = sin(phi);
103  a[j] = EGS_Vector(-sphi,cphi,0)*R;
104  d[j] = a[j]*xo;
105  }
106 }
107 
108 EGS_IPlanes::~EGS_IPlanes() {
109  delete [] a;
110  delete [] d;
111 }
112 
113 int EGS_IPlanes::isWhere(const EGS_Vector &x) {
114  EGS_Float aux_old = a[0]*x;
115  for (int j=1; j<nreg; j++) {
116  EGS_Float aux = a[j]*x;
117  if (aux_old >= d[j-1] && aux < d[j]) {
118  return j-1;
119  }
120  aux_old = aux;
121  }
122  return nreg-1;
123 }
124 
125 int EGS_IPlanes::inside(const EGS_Vector &x) {
126  EGS_Float aux_old = a[0]*x;
127  for (int j=1; j<nreg; j++) {
128  EGS_Float aux = a[j]*x;
129  if (aux_old >= d[j-1] && aux < d[j]) {
130  return j-1;
131  }
132  aux_old = aux;
133  }
134  return nreg-1;
135 }
136 
137 int EGS_IPlanes::howfar(int ireg, const EGS_Vector &x,
138  const EGS_Vector &u, EGS_Float &t, int *newmed, EGS_Vector *normal) {
139  if (ireg < 0) {
140  egsFatal("EGS_IPlanes::howfar: ireg (%d) can not be negative\n",ireg);
141  }
142  EGS_Float t1 = t, t2 = t;
143  EGS_Float up = a[ireg]*u;
144  int inew = ireg;
145  if (up < 0) {
146  t1 = (d[ireg] - a[ireg]*x)/up;
147  if (t1 <= t) {
148  t = t1;
149  inew = ireg-1;
150  if (inew < 0) {
151  inew = nreg-1;
152  }
153  if (newmed) {
154  *newmed = medium(inew);
155  }
156  if (normal) {
157  *normal = a[ireg];
158  }
159  }
160  }
161  int j = ireg+1;
162  if (j >= nreg) {
163  j = 0;
164  }
165  up = a[j]*u;
166  if (up > 0) {
167  t2 = (d[j] - a[j]*x)/up;
168  if (t2 < t) {
169  t = t2;
170  inew = j;
171  if (newmed) {
172  *newmed = medium(inew);
173  }
174  if (normal) {
175  *normal = a[j]*(-1);
176  }
177  }
178  }
179  return inew;
180 }
181 
182 EGS_Float EGS_IPlanes::hownear(int ireg, const EGS_Vector &x) {
183  if (ireg < 0) {
184  egsFatal("EGS_IPlanes::hownear: ireg (%d) can not be negative\n",ireg);
185  }
186  EGS_Float t1 = a[ireg]*x - d[ireg];
187  int j = ireg+1;
188  if (j >= nreg) {
189  j = 0;
190  }
191  EGS_Float t2 = d[j] - a[j]*x;
192  return min(t1,t2);
193 }
194 
195 void EGS_IPlanes::printInfo() const {
197  egsInformation(" axis: Xo = (%g,%g,%g) a = (%g,%g,%g)\n",xo.x,xo.y,xo.z,
198  axis.x,axis.y,axis.z);
200  "\n=======================================================\n");
201 }
202 
203 /***************************************************************************/
204 
205 EGS_RadialRepeater::EGS_RadialRepeater(const EGS_Vector &Xo,
206  const EGS_Vector &A, int np, EGS_BaseGeometry *G,
207  EGS_Float first, const string &Name) : EGS_BaseGeometry(Name), g(G) {
208  EGS_Float dphi = 2*M_PI/np;
209  iplanes = new EGS_IPlanes(Xo,A,np,first-dphi/2);
210  //iplanes = new EGS_IPlanes(Xo,A,np,first);
211  //iplanes = new EGS_IPlanes(Xo,A,np,(EGS_Float)0);
212  g->ref();
213  iplanes->ref();
214  nrep = np;
215  ng = g->regions();
216  nreg = nrep*ng + 1;
217  //EGS_Float dphi = 2*M_PI/np;
218  R = new EGS_RotationMatrix [nrep];
219  EGS_RotationMatrix Ro(A);
220  for (int j=0; j<nrep; j++) {
221  //EGS_Float phi = first + dphi*(0.5+j);
222  //EGS_Float phi = dphi*(0.5+j);
223  EGS_Float phi = dphi*j;
224  R[j] = Ro.inverse()*EGS_RotationMatrix::rotZ(-phi)*Ro;
225  }
226  phi_o = first;
227  xo = Xo;
228 }
229 
230 EGS_RadialRepeater::~EGS_RadialRepeater() {
231  delete [] R;
232  if (!iplanes->deref()) {
233  delete iplanes;
234  }
235  if (!g->deref()) {
236  delete g;
237  }
238 }
239 
240 void EGS_RadialRepeater::printInfo() const {
242  egsInformation("%d uniformely rotated replicas of geometry %s with "
243  "phi_o=%g degrees\n",nrep,g->getName().c_str(),
244  phi_o*180./M_PI);
245  EGS_Vector xo(iplanes->getAxisXo()), axis(iplanes->getAxisDirection());
246  egsInformation("Axis of rotation is xo=(%g,%g,%g) a=(%g,%g,%g)\n",
247  xo.x,xo.y,xo.z,axis.x,axis.y,axis.z);
248  const char *med_name = getMediumName(med);
249  if (med_name) {
250  egsInformation("space is filled with %s\n",med_name);
251  }
252  else {
253  egsInformation("space is filled with vacuum\n");
254  }
255  egsInformation("Repeated geometry:\n");
256  g->printInfo();
257 }
258 
259 extern "C" {
260 
261  EGS_IPLANES_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
262  if (!input) {
263  egsWarning("createGeometry(iplanes): null input?\n");
264  return 0;
265  }
266  vector<EGS_Float> axis;
267  EGS_Vector xo, a(0,0,1);
268  int err = input->getInput("axis",axis);
269  if (!err && axis.size() == 6) {
270  xo = EGS_Vector(axis[0],axis[1],axis[2]);
271  a = EGS_Vector(axis[3],axis[4],axis[5]);
272  a.normalize();
273  }
274  else egsWarning("createGeometry(iplanes): wrong/missing axis input\n"
275  " using Xo=(0,0,0), a=(0,0,1)\n");
276  string type;
277  err = input->getInput("type",type);
278  if (!err && (type == "EGS_RadialRepeater" || type == "repeater")) {
279  int np;
280  err = input->getInput("number of repetitions",np);
281  if (err || np < 2) {
282  egsWarning("createGeometry(iplanes): missing/wrong "
283  "'number of repetitions' input\n");
284  return 0;
285  }
286  string gname;
287  err = input->getInput("repeated geometry",gname);
288  if (err) {
289  egsWarning("createGeometry(iplanes): missing 'repeated geometry'"
290  " input\n");
291  return 0;
292  }
294  if (!g) {
295  egsWarning("createGeometry(iplanes): no geometry named %s exists\n",
296  gname.c_str());
297  return 0;
298  }
299  EGS_Float phi_o = 0;
300  EGS_Float tmp1,tmp2;
301  err = input->getInput("first angle",tmp1);
302  int err1 = input->getInput("first angle in radians",tmp2);
303  if (!err) {
304  phi_o = tmp1*M_PI/180;
305  }
306  else if (!err1) {
307  phi_o = tmp2;
308  }
309  EGS_RadialRepeater *result = new EGS_RadialRepeater(xo,a,np,g,phi_o);
310  result->setName(input);
311  result->setBoundaryTolerance(input);
312  string medium;
313  err = input->getInput("medium",medium);
314  if (!err) {
315  result->setMedium(medium);
316  }
317  result->setRLabels(input);
318  result->setLabels(input);
319  return result;
320  }
321 
322  vector<EGS_Float> angles;
323  err = input->getInput("angles",angles);
324  bool is_degree = true;
325  EGS_Float max_angle = 180;
326  if (err) {
327  err = input->getInput("angles in radian",angles);
328  if (err) {
329  egsWarning("createGeometry(iplanes): wrong/missing 'angles' or "
330  "'angles in radian' input\n");
331  return 0;
332  }
333  is_degree = false;
334  max_angle = M_PI;
335  }
336  EGS_Float *ang = new EGS_Float [angles.size()];
337  for (int j=0; j<angles.size(); j++) {
338  ang[j] = angles[j];
339  /*
340  if( ang[j] < 0 || ang[j] > max_angle ) {
341  egsWarning("createGeometry(iplanes): angle must be between 0 and "
342  "%g\n",max_angle); delete [] ang; return 0;
343  }
344  */
345  if (j > 0) {
346  if (ang[j] <= ang[j-1]) {
347  egsWarning("createGeometry(iplanes): angles must be ordered"
348  " in increasing order\n");
349  delete [] ang;
350  return 0;
351  }
352  }
353  }
354  int nang = angles.size();
355  EGS_Float ang_diff = ang[nang-1]-ang[0];
356  if (ang_diff > max_angle) {
357  egsWarning("createGeometry(iplanes): difference between first and last"
358  " angle must be less then 180 degrees.\n");
359  if (is_degree)
360  egsWarning(" Your input: first angle=%g, last angle=%g,"
361  " difference=%g\n",ang[0],ang[nang-1],ang[nang-1]-ang[0]);
362  else
363  egsWarning(" Your input: first angle=%g, last angle=%g,"
364  " difference=%g\n",ang[0]*180/M_PI,ang[nang-1]*180/M_PI,
365  (ang[nang-1]-ang[0])*180/M_PI);
366  delete [] ang;
367  return 0;
368  }
369  EGS_BaseGeometry *g = new EGS_IPlanes(xo,a,angles.size(),ang,"",is_degree);
370  delete [] ang;
371  g->setName(input);
372  g->setBoundaryTolerance(input);
373  g->setMedia(input);
374  g->setLabels(input);
375  return g;
376  }
377 
378  void EGS_RadialRepeater::setRLabels(EGS_Input *input) {
379 
380  // radial repeater labels
381  string inp;
382  int err;
383  err = input->getInput("set repeater label", inp);
384  if (!err) {
385  iplanes->setLabels(inp);
386  }
387  }
388 
389 
390  void EGS_RadialRepeater::getLabelRegions(const string &str, vector<int> &regs) {
391 
392  vector<int> local_regs;
393 
394  // label in repeated geometry
395  local_regs.clear();
396  g->getLabelRegions(str, local_regs);
397  for (int i=0; i<nrep; i++) {
398  for (int r=0; r<local_regs.size(); r++) {
399  regs.push_back(ng*i + local_regs[r]);
400  }
401  }
402 
403  // labels defined in iplanes
404  local_regs.clear();
405  iplanes->getLabelRegions(str, local_regs);
406  for (int i=0; i<local_regs.size(); i++) {
407  for (int r=0; r<ng; r++) {
408  regs.push_back(ng*local_regs[i] + r);
409  }
410  }
411 
412  // label defined in self (repeater input block)
414 
415  }
416 
417 }
int deref()
Decrease the reference count to this geometry.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
EGS_IPlanes(const EGS_Vector &Xo, const EGS_Vector &A, int np, const EGS_Float *angles, const string &Name="", bool degree=true)
Construct a set of intersecting planes (iplanes)
Definition: egs_iplanes.cpp:50
A class for vector rotations.
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
int setLabels(EGS_Input *input)
Set the labels from an input block.
A set of planes intersecting in the same axis.
Definition: egs_iplanes.h:150
Global egspp functions header file.
int med
Medium index.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
EGS_Float z
z-component
Definition: egs_vector.h:62
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
static EGS_RotationMatrix rotZ(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the z-axis by the angle with cphi, sphi = .
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
void setMedium(const string &Name)
Set all regions to a medium with name Name.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
EGS_AffineTransform and EGS_RotationMatrix class header file.
Attempts to fix broken math header files.
A radial geometry replicator.
Definition: egs_iplanes.h:341
EGS_Float x
x-component
Definition: egs_vector.h:60
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
static const char * getMediumName(int ind)
Get the name of medium with index ind.
const string & getName() const
Get the name of this geometry.
Intersecting planes: header.
int nreg
Number of local regions in this geometry.
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.