EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_cones.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ cones 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_cones.h"
38 #include "egs_input.h"
39 
40 #include <vector>
41 using std::vector;
42 
43 #ifndef M_PI
44  #define M_PI 3.14159265358979323846
45 #endif
46 
47 string EGS_SimpleCone::type = "EGS_SimpleCone";
48 string EGS_ParallelCones::type = "EGS_ParallelCones";
49 string EGS_ConeSet::type = "EGS_ConeSet";
50 string EGS_ConeStack::type = "EGS_ConeStack";
51 
52 void EGS_ConeStack::clear(bool all) {
53  if (nltot > 0) {
54  if (all) {
55  for (int j=0; j<nl; j++) {
56  for (int i=0; i<nr[j]; i++)
57  if (!cones[j][i]->deref()) {
58  delete cones[j][i];
59  }
60  delete [] cones[j];
61  }
62  nl = 0;
63  }
64  delete [] cones;
65  delete [] pos;
66  delete [] nr;
67  delete [] flag;
68  nltot = 0;
69  }
70 }
71 
72 #define N_CS_GROW 50
73 void EGS_ConeStack::resize() {
74  int nnew = nltot + N_CS_GROW;
75  int *new_nr = new int [nnew], *new_flag = new int [nnew];
76  EGS_Float *new_pos = new EGS_Float [nnew+1];
77  EGS_SimpleCone ***new_cones = new EGS_SimpleCone **[nnew];
78  for (int j=0; j<nl; j++) {
79  new_nr[j] = nr[j];
80  new_flag[j] = flag[j];
81  new_pos[j] = pos[j];
82  new_cones[j] = cones[j];
83  }
84  if (nl > 0) {
85  new_pos[nl] = pos[nl];
86  }
87  clear(false);
88  nr = new_nr;
89  pos = new_pos;
90  flag = new_flag;
91  cones = new_cones;
92  nltot = nnew;
93 }
94 
95 void EGS_ConeStack::addLayer(EGS_Float thick, const vector<EGS_Float> &rtop,
96  const vector<EGS_Float> &rbottom,
97  const vector<string> &med_names) {
98  int this_nr = rbottom.size();
99  if (!this_nr) {
100  egsWarning("EGS_ConeStack::addLayer: no bottom radii?\n");
101  egsWarning(" --> ignoring layer\n");
102  return;
103  }
104  if ((int)med_names.size() != this_nr) {
105  egsWarning("EGS_ConeStack::addLayer: number of cone radii (%d) is"
106  " different from number of media (%d)\n",this_nr,med_names.size());
107  egsWarning(" --> ignoting layer\n");
108  return;
109  }
110  bool same_radii = false;
111  if (rtop.size() != rbottom.size()) {
112  if (!rtop.size()) {
113  if (!nl) {
114  egsWarning("EGS_ConeStack::addLayer: zero top radii does not"
115  " work for the first layer\n");
116  egsWarning(" --> ignoting layer\n");
117  return;
118  }
119  same_radii = true;
120  }
121  else {
122  egsWarning("EGS_ConeStack::addLayer: number of bottom radii (%d)"
123  " is different from number of top radii (%d)\n",this_nr,
124  rtop.size());
125  egsWarning(" --> ignoting layer\n");
126  return;
127  }
128  }
129  if (nl >= nltot) {
130  resize();
131  }
132  if (!nl) {
133  pos[nl] = xo*a;
134  Rout = rbottom[this_nr-1];
135  Rout2 = Rout*Rout;
136  if (fabs(rtop[this_nr-1]-Rout) > boundaryTolerance) {
137  same_Rout = false;
138  is_convex = false;
139  }
140  }
141 
142  cones[nl] = new EGS_SimpleCone * [this_nr];
143  nr[nl] = this_nr;
144  pos[nl+1] = pos[nl] + thick;
145  EGS_Vector x(xo+a*(pos[nl]-pos[0]));
146  //egsInformation(" layer %d x = (%g,%g,%g)\n",nl,x.x,x.y,x.z);
147  for (int ir=0; ir<this_nr; ir++) {
148  EGS_Float Rtop = same_radii ? cones[nl-1][ir]->getRadius(x) : rtop[ir];
149  cones[nl][ir] = new EGS_SimpleCone(x,a,thick,Rtop,rbottom[ir]);
150  cones[nl][ir]->setMedium(med_names[ir]);
151  cones[nl][ir]->ref();
152  if (ir == this_nr-1) {
153  if (fabs(rbottom[this_nr-1]-Rout) > boundaryTolerance) {
154  same_Rout = false;
155  is_convex = false;
156  }
157  if (fabs(Rtop-Rout) > boundaryTolerance) {
158  same_Rout = false;
159  is_convex = false;
160  }
161  }
162  }
163  if (same_radii) {
164  flag[nl-1] += 2;
165  flag[nl] = 1;
166  }
167  else {
168  flag[nl] = 0;
169  }
170  if (this_nr > nmax) {
171  nmax = this_nr;
172  }
173  ++nl;
174  nreg = nl*nmax;
175  //egsInformation("addLayer: the following layers are defined:\n");
176  //for(int j=0; j<nl; j++) egsInformation(" %g %g %d %d\n",pos[j],pos[j+1],
177  // nr[j],flag[j]);
178 }
179 
180 void EGS_ConeStack::printInfo() const {
182  egsInformation("number of layers: %d\n",nl);
183  for (int il=0; il<nl; il++) {
184  egsInformation("*** layer %d: top = %g bottom = %g\n",il,
185  pos[il],pos[il+1]);
186  egsInformation(" top radii: ");
187  int i;
188  EGS_Vector x(xo+a*(pos[il]-pos[0]));
189  for (i=0; i<nr[il]; i++) {
190  egsInformation("%g ",cones[il][i]->getRadius(x));
191  }
192  egsInformation("\n bottom radii: ");
193  x = xo + a*(pos[il+1]-pos[0]);
194  for (i=0; i<nr[il]; i++) {
195  egsInformation("%g ",cones[il][i]->getRadius(x));
196  }
197  egsInformation("\n media: ");
198  for (i=0; i<nr[il]; i++) {
199  egsInformation("%d ",cones[il][i]->medium(0));
200  }
201  egsInformation("\n");
202  }
203  egsInformation("=====================================================\n");
204 }
205 
206 void EGS_SimpleCone::printInfo() const {
208  egsInformation(" cone apex = (%g,%g,%g)\n",xo.x,xo.y,xo.z);
209  egsInformation(" cone axis = (%g,%g,%g)\n",a.x,a.y,a.z);
210  egsInformation(" opening angle = %g degrees\n",180*atan(gamma)/M_PI);
211  if (open) {
212  egsInformation(" cone is open\n");
213  }
214  else {
215  egsInformation(" cone height %g\n",-d1);
216  }
217  egsInformation("=====================================================\n");
218 }
219 
220 void EGS_ParallelCones::printInfo() const {
222  egsInformation(" cone axis = (%g,%g,%g)\n",a.x,a.y,a.z);
223  egsInformation(" opening angle = %g degrees\n",180*atan(gamma)/M_PI);
224  egsInformation(" number of cones = %d\n",nc);
225  egsInformation(" cone apexes:\n ");
226  for (int j=0; j<nc; j++) {
227  egsInformation("(%g,%g,%g) ",xo[j].x,xo[j].y,xo[j].z);
228  }
229  egsInformation("\n=====================================================\n");
230 }
231 
232 void EGS_ConeSet::printInfo() const {
234  egsInformation(" cone apex = (%g,%g,%g)\n",xo.x,xo.y,xo.z);
235  egsInformation(" cone axis = (%g,%g,%g)\n",a.x,a.y,a.z);
236  egsInformation(" flag = %d\n",flag);
237  egsInformation(" opening angles (degrees)=\n ");
238  for (int j=0; j<nc; j++) {
239  egsInformation("%g ",180*atan(gamma[j])/M_PI);
240  }
241  egsInformation("\n=====================================================\n");
242 }
243 
244 
245 extern "C" {
246 
247  EGS_CONES_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
248 
249  if (!input) {
250  egsWarning("createGeometry(cones): null input?\n");
251  return 0;
252  }
253  string type;
254  int err = input->getInput("type",type);
255  if (err) {
256  egsWarning("createGeometry(cones): object type not specified\n");
257  return 0;
258  }
259  if (type == "EGS_ConeStack") {
260  vector<EGS_Float> axis;
261  err = input->getInput("axis",axis);
262  if (err) {
263  egsWarning("createGeometry(EGS_ConeStack): no axis input\n");
264  return 0;
265  }
266  if (axis.size() != 6) {
267  egsWarning("createGeometry(EGS_ConeStack): wrong axis input"
268  " (expecting 6 instead of %d inputs)\n",axis.size());
269  return 0;
270  }
271  EGS_ConeStack *g = new EGS_ConeStack(
272  EGS_Vector(axis[0],axis[1],axis[2]),
273  EGS_Vector(axis[3],axis[4],axis[5]),"");
274  EGS_Input *layer;
275  int nl = 0;
276  vector<int> layerLabels;
277  while ((layer = input->takeInputItem("layer"))) {
278  vector<EGS_Float> rtop, rbottom;
279  vector<string> media;
280  EGS_Float thickness;
281  err = layer->getInput("thickness",thickness);
282  if (err)
283  egsWarning("createGeometry(EGS_ConeStack): missing 'thickness'"
284  " input for layer %d\n --> layer ignored\n",nl);
285  else {
286  err = layer->getInput("top radii",rtop);
287  err = layer->getInput("bottom radii",rbottom);
288  if (err) egsWarning("createGeometry(EGS_ConeStack): missing "
289  "'bottom radii' input for layer %d\n",nl);
290  int err1 = layer->getInput("media",media);
291  if (err1) egsWarning("createGeometry(EGS_ConeStack): missing "
292  "'media' input for layer %d\n",nl);
293  if (err || err1) {
294  egsWarning(" --> layer ignored\n");
295  }
296  else {
297  g->addLayer(thickness,rtop,rbottom,media);
298  }
299  }
300  layerLabels.push_back(g->setLabels(layer));
301  delete layer;
302  ++nl;
303  }
304  if (!g->nLayer()) {
305  egsWarning("createGeometry(EGS_ConeStack): zero layers\n");
306  delete g;
307  return 0;
308  }
309 
310  // adjust lable region numbering in each layer
311  int count=0;
312  for (size_t i=0; i<layerLabels.size(); i++) {
313  for (int j=0; j<layerLabels[i]; j++) {
314  g->shiftLabelRegions(count,i);
315  count++;
316  }
317  }
318 
319  g->setName(input);
320  g->setBoundaryTolerance(input);
321  g->setLabels(input);
322  g->setBScaling(input); // Perhaps add density scaling as well?
323  return g;
324  }
325 
326  // get cone apex position
327  vector<EGS_Float> tmp;
328  EGS_Vector Xo;
329  err = input->getInput("apex",tmp);
330  if (err || tmp.size() != 3)
331  egsWarning("createGeometry(cones): no 'apex' input, "
332  "assuming (0,0,0)\n");
333  else {
334  Xo.x = tmp[0];
335  Xo.y = tmp[1];
336  Xo.z = tmp[2];
337  }
338  tmp.clear();
339 
340  // get cone axis
341  EGS_Vector axis(0,0,1);
342  err = input->getInput("axis",tmp);
343  if (err || tmp.size() != 3)
344  egsWarning("createGeometry(cones): no 'axis' input, assuming "
345  "(0,0,1)\n");
346  else {
347  axis.x = tmp[0];
348  axis.y = tmp[1];
349  axis.z = tmp[2];
350  }
351 
352  EGS_BaseGeometry *g;
353  if (input->compare(type,"EGS_ConeSet")) {
354  vector<EGS_Float> angles;
355  err = input->getInput("opening angles",angles);
356  bool is_radian = false;
357  if (err) {
358  angles.clear();
359  err = input->getInput("opening angles in radian",angles);
360  if (err) {
361  egsWarning("createGeometry(cones): no 'opening angles' or "
362  "'opening angles in radian' input\n");
363  return 0;
364  }
365  is_radian = true;
366  }
367  int flag = 0;
368  err = input->getInput("flag",flag);
369  int nc = angles.size();
370  EGS_Float *gamma = new EGS_Float [nc];
371  for (int j=0; j<nc; j++) {
372  if (angles[j] <= 0) {
373  egsWarning("createGeometry(cones): opening angles must be"
374  " positive\n");
375  delete [] gamma;
376  return 0;
377  }
378  EGS_Float a = angles[j];
379  if (!is_radian) {
380  a *= (M_PI/180);
381  }
382  if (a >= M_PI/2) {
383  egsWarning("createGeometry(cones): opening angles can not be"
384  " greater than Pi/2\n");
385  delete [] gamma;
386  return 0;
387  }
388  if (j > 0) {
389  if (angles[j] <= angles[j-1]) {
390  egsWarning("createGeometry(cones): opening angles must be"
391  " in increasing order\n");
392  delete [] gamma;
393  return 0;
394  }
395  }
396  gamma[j] = tan(a);
397  }
398  g = new EGS_ConeSet(Xo,axis,nc,gamma,flag);
399  delete [] gamma;
400  }
401 
402  else {
403  // get opening angle
404  EGS_Float angle;
405  err = input->getInput("opening angle",angle);
406  if (err) {
407  err = input->getInput("opening angle in radian",angle);
408  if (err) {
409  egsWarning("createGeometry(cones): no 'opening angle' or "
410  "'opening angle in radian' input\n");
411  return 0;
412  }
413  }
414  else {
415  angle *= (M_PI/180);
416  }
417  if (angle >= M_PI/2) {
418  egsWarning("createGeometry(cones): it is not allowed to have"
419  " an opening angle greater than Pi/2, your input was %g\n",angle);
420  return 0;
421  }
422  EGS_Float gamma = tan(angle);
423  if (input->compare(type,"EGS_SimpleCone")) {
424  EGS_Float height;
425  EGS_Float *h = 0;
426  err = input->getInput("height",height);
427  if (!err) {
428  if (height <= 0)
429  egsWarning("createGeometry(cones): the cone height"
430  " must be greater than zero, your input was %g\n",height);
431  else {
432  h = &height;
433  }
434  }
435  g = new EGS_SimpleCone(Xo,axis,gamma,h,"");
436  }
437  else if (input->compare(type,"EGS_ParallelCones")) {
438  vector<EGS_Float> d;
439  err = input->getInput("apex distances",d);
440  EGS_Float *dist=0;
441  int nc=1;
442  if (!err && d.size() > 0) {
443  dist = new EGS_Float [d.size()];
444  for (size_t j=0; j<d.size(); j++) {
445  dist[j] = d[j];
446  }
447  nc = 1 + d.size();
448  }
449  g = new EGS_ParallelCones(nc,Xo,axis,gamma,dist);
450  if (dist) {
451  delete [] dist;
452  }
453  }
454  else {
455  egsWarning("createGeometry(cones): unknown object type %s\n",
456  type.c_str());
457  return 0;
458  }
459  }
460 
461  g->setName(input);
462  g->setMedia(input);
463  g->setLabels(input);
464  return g;
465  }
466 
467 }
A single cone that may be open (i.e. extends to infinity or closed by a plane perpendicular to the co...
Definition: egs_cones.h:124
int deref()
Decrease the reference count to this geometry.
virtual void printInfo() const
Print information about this geometry.
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 &quot;parallel cones&quot; (i.e. cones with the same axis and opening angles but different apexes) ...
Definition: egs_cones.h:739
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.
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
EGS_Float z
z-component
Definition: egs_vector.h:62
Various cone geometries: header.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
A set of cones with different opening angles but the same axis and apexes.
Definition: egs_cones.h:1059
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
bool is_convex
Is this geometry convex?
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_Float boundaryTolerance
Boundary tolerance for geometries that need it.
EGS_Float x
x-component
Definition: egs_vector.h:60
A cone stack.
Definition: egs_cones.h:1461
int ref()
Increase the reference count to this geometry.
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 bool compare(const string &s1, const string &s2)
Definition: egs_input.cpp:1170
int nreg
Number of local regions in this geometry.
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
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.