EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_envelope_geometry.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ envelope 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: Ernesto Mainegra-Hing
27 # Frederic Tessier
28 # Reid Townson
29 # Max Orok
30 #
31 ###############################################################################
32 */
33 
34 
40 #ifndef EGS_ENVELOPE_GEOMETRY_
41 #define EGS_ENVELOPE_GEOMETRY_
42 
43 #ifdef WIN32
44 
45  #ifdef BUILD_ENVELOPEG_DLL
46  #define EGS_ENVELOPEG_EXPORT __declspec(dllexport)
47  #else
48  #define EGS_ENVELOPEG_EXPORT __declspec(dllimport)
49  #endif
50  #define EGS_ENVELOPEG_LOCAL
51 
52 #else
53 
54  #ifdef HAVE_VISIBILITY
55  #define EGS_ENVELOPEG_EXPORT __attribute__ ((visibility ("default")))
56  #define EGS_ENVELOPEG_LOCAL __attribute__ ((visibility ("hidden")))
57  #else
58  #define EGS_ENVELOPEG_EXPORT
59  #define EGS_ENVELOPEG_LOCAL
60  #endif
61 
62 #endif
63 
64 
65 #include "egs_base_geometry.h"
66 #include "egs_functions.h"
67 
68 #include<vector>
69 using std::vector;
70 
229 class EGS_ENVELOPEG_EXPORT EGS_EnvelopeGeometry : public EGS_BaseGeometry {
230 
231 public:
232 
234  const vector<EGS_BaseGeometry *> &geoms, const string &Name = "",
235  bool newindexing=false);
236 
238 
239  bool isRealRegion(int ireg) const {
240  if (ireg < 0 || ireg >= nreg) {
241  return false;
242  }
243  if (ireg < nbase) {
244  return g->isRealRegion(ireg);
245  }
246  int jg, ilocal;
247  if (new_indexing) {
248  jg = reg_to_inscr[ireg-nbase];
249  ilocal = ireg - local_start[jg];
250  }
251  else {
252  jg = (ireg - nbase)/nmax;
253  ilocal = ireg - nbase - jg*nmax;
254  }
255  return geometries[jg]->isRealRegion(ilocal);
256  };
257 
258  bool isInside(const EGS_Vector &x) {
259  return g->isInside(x);
260  };
261 
262  int isWhere(const EGS_Vector &x) {
263  int ireg = g->isWhere(x);
264  if (ireg < 0) {
265  return ireg;
266  }
267  for (int j=0; j<n_in; j++) {
268  int i = geometries[j]->isWhere(x);
269  if (i >= 0) return new_indexing ? local_start[j] + i :
270  nbase + nmax*j + i;
271  }
272  return ireg;
273  };
274 
275  int inside(const EGS_Vector &x) {
276  return isWhere(x);
277  };
278 
279  int medium(int ireg) const {
280  if (ireg < nbase) {
281  return g->medium(ireg);
282  }
283  int jg, ilocal;
284  if (new_indexing) {
285  jg = reg_to_inscr[ireg-nbase];
286  ilocal = ireg - local_start[jg];
287  }
288  else {
289  jg = (ireg - nbase)/nmax;
290  ilocal = ireg - nbase - jg*nmax;
291  }
292  return geometries[jg]->medium(ilocal);
293  };
294 
295  int computeIntersections(int ireg, int n, const EGS_Vector &X,
296  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
297  if (n < 1) {
298  return -1;
299  }
300  int ifirst = 0;
301  EGS_Float t, ttot = 0;
302  EGS_Vector x(X);
303  int imed;
304  if (ireg < 0) {
305  t = veryFar;
306  ireg = howfar(ireg,x,u,t,&imed);
307  if (ireg < 0) {
308  return 0;
309  }
310  isections[0].t = t;
311  isections[0].rhof = 1;
312  isections[0].ireg = -1;
313  isections[0].imed = -1;
314  ttot = t;
315  ++ifirst;
316  x += u*t;
317  }
318  else {
319  imed = medium(ireg);
320  }
321 
322 
323  int j = ifirst;
324  int ij = -1, ig=0;
325  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
326  if (loopCount == loopMax) {
327  egsFatal("EGS_EnvelopeGeometry::computeIntersections: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
328  return -1;
329  }
330  if (ireg >= nbase) {
331  if (new_indexing) {
332  ig = reg_to_inscr[ireg-nbase];
333  ij = ireg - local_start[ig];
334  }
335  else {
336  ig = (ireg - nbase)/nmax;
337  ij = ireg - nbase - ig*nmax;
338  }
339  }
340  isections[j].imed = imed;
341  isections[j].ireg = ireg;
342  isections[j].rhof = getRelativeRho(ireg);
343  if (ireg < nbase) { // in one of the regions of the base geometry
344  t = veryFar;
345  int ibase = g->howfar(ireg,x,u,t,&imed);
346  ij = -1;
347  for (int i=0; i<n_in; i++) {
348  int ireg_i = geometries[i]->howfar(-1,x,u,t,&imed);
349  if (ireg_i >= 0) {
350  ij = ireg_i;
351  ig = i;
352  }
353  }
354  ttot += t;
355  isections[j++].t = ttot;
356  if (ij < 0) {
357  ireg = ibase;
358  }
359  else ireg = new_indexing ? local_start[ig] + ij :
360  nbase + ig*nmax + ij;
361  if (ireg < 0) {
362  return j;
363  }
364  if (j >= n) {
365  return -1;
366  }
367  x += u*t;
368  }
369  else {
370  int iadd = new_indexing ? local_start[ig] : nbase + ig*nmax;
371  int nsec = geometries[ig]->computeIntersections(ij,n-j,
372  x,u,&isections[j]);
373  int nm = nsec >= 0 ? nsec+j : n;
374  for (int i=j; i<nm; i++) {
375  isections[i].ireg += iadd;
376  isections[i].t += ttot;
377  }
378  if (nsec < 0) {
379  return nsec;
380  }
381  j += nsec;
382  if (j >= n) {
383  return -1;
384  }
385  t = isections[j-1].t - ttot;
386  x += u*t;
387  ttot = isections[j-1].t;
388  ireg = g->isWhere(x);
389  if (ireg < 0) {
390  return j;
391  }
392  imed = g->medium(ireg);
393  }
394  }
395  return -1;
396  }
397 
398  EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
399  const EGS_Vector &u) {
400  if (ireg < 0) {
401  return 0;
402  }
403  EGS_Float d;
404  if (ireg < nbase) {
405  d = g->howfarToOutside(ireg,x,u);
406  }
407  else if (g->regions() == 1) {
408  d = g->howfarToOutside(0,x,u);
409  }
410  else {
411  int ir = g->isWhere(x);
412  d = g->howfarToOutside(ir,x,u);
413  }
414  return d;
415  };
416 
417  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
418  EGS_Float &t, int *newmed = 0, EGS_Vector *normal = 0) {
419  if (ireg >= 0) {
420  // inside.
421  if (ireg < nbase) {
422  // in one of the regions of the base geometry
423  // check if we hit a boundary in the base geometry.
424  // if we do, newmed and normal get set accordingly.
425  int ibase = g->howfar(ireg,x,u,t,newmed,normal);
426  int ij = -1, jg;
427  // check if we will enter any of the inscribed geometries
428  // before entering a new region in the base geometry.
429  for (int j=0; j<n_in; j++) {
430  int ireg_j =
431  geometries[j]->howfar(-1,x,u,t,newmed,normal);
432  if (ireg_j >= 0) {
433  // we do. remember the inscribed geometry index
434  // and local region
435  ij = ireg_j;
436  jg = j;
437  }
438  }
439  if (ij < 0) {
440  return ibase;
441  }
442  // ij<0 implies that we have not hit any of the
443  // inscribed geometries => return the base geometry index.
444  // ij>=0 implies that we entered inscribed geometry
445  // jg in its local region ij.
446  return new_indexing ? local_start[jg] + ij :
447  nbase + jg*nmax + ij;
448  }
449  // if here, we are in an inscribed geometry.
450  // calculate its index (jg) and its local region (ilocal).
451  int jg, ilocal;
452  if (new_indexing) {
453  jg = reg_to_inscr[ireg-nbase];
454  ilocal = ireg-local_start[jg];
455  }
456  else {
457  jg = (ireg - nbase)/nmax;
458  ilocal = ireg - nbase - jg*nmax;
459  }
460  // and then check if we will hit a boundary in this geometry.
461  int inew = geometries[jg]->howfar(ilocal,x,u,t,newmed,normal);
462  if (inew >= 0) return new_indexing ? local_start[jg] + inew :
463  nbase + jg*nmax + inew;
464  // inew >= 0 implies that we either stay in the same
465  // region (inew=ilocal) or we entered a new region
466  // (inew!=ilocal), which is still inside the inscribed geometry
467  // inew<0 implies that we have exited the inscribed geometry
468  // => check to see in which base geometry region we are.
469  inew = g->isWhere(x+u*t);
470  if (inew >= 0 && newmed) {
471  *newmed = g->medium(inew);
472  }
473  return inew;
474  }
475  // if here, we are outside the base geometry.
476  // check to see if we will enter.
477  int ienter = g->howfar(ireg,x,u,t,newmed,normal);
478  if (ienter >= 0) {
479  // yes, we do. see if we are already inside of one of the
480  // inscribed geometries.
481  for (int j=0; j<n_in; j++) {
482  int i = geometries[j]->isWhere(x+u*t);
483  if (i >= 0) {
484  // yes, we are.
485  if (newmed) {
486  *newmed = geometries[j]->medium(i);
487  }
488  return new_indexing ? local_start[j] + i :
489  nbase + nmax*j + i;
490  }
491  }
492  }
493  return ienter;
494  };
495 
496  EGS_Float hownear(int ireg, const EGS_Vector &x) {
497  if (ireg >= 0) {
498  EGS_Float tmin;
499  if (ireg < nbase) { // in one of the regions of the base geom.
500  tmin = g->hownear(ireg,x);
501  for (int j=0; j<n_in; j++) {
502  EGS_Float tj = geometries[j]->hownear(-1,x);
503  if (tj < tmin) {
504  tmin = tj;
505  if (tmin <= 0) {
506  return tmin;
507  }
508  }
509  }
510  return tmin;
511  }
512  int jg, ilocal;
513  if (new_indexing) {
514  jg = reg_to_inscr[ireg-nbase];
515  ilocal = ireg-local_start[jg];
516  }
517  else {
518  jg = (ireg - nbase)/nmax;
519  ilocal = ireg - nbase - jg*nmax;
520  }
521  return geometries[jg]->hownear(ilocal,x);
522  }
523  return g->hownear(ireg,x);
524  };
525 
526  int getMaxStep() const {
527  int nstep = g->getMaxStep();
528  for (int j=0; j<n_in; ++j) {
529  nstep += geometries[j]->getMaxStep();
530  }
531  return nstep+n_in;
532  };
533 
534  bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
535  if (ireg >= 0 && ireg < nreg) {
536  if (ireg < nbase) {
537  return g->hasBooleanProperty(ireg,prop);
538  }
539  int jg = (ireg - nbase)/nmax;
540  int ilocal = ireg - nbase - jg*nmax;
541  return geometries[jg]->hasBooleanProperty(ilocal,prop);
542  }
543  return false;
544  };
545  void setBooleanProperty(EGS_BPType) {
546  setPropertyError("setBooleanProperty()");
547  };
548  void addBooleanProperty(int) {
549  setPropertyError("addBooleanProperty()");
550  };
551  void setBooleanProperty(EGS_BPType,int,int,int step=1) {
552  setPropertyError("setBooleanProperty()");
553  };
554  void addBooleanProperty(int,int,int,int step=1) {
555  setPropertyError("addBooleanProperty()");
556  };
557 
558  const string &getType() const {
559  return type;
560  };
561 
562  EGS_BaseGeometry **getInscribedGeometries(std::size_t &nInscribed) const {
563  nInscribed = n_in;
564  return geometries;
565  }
566 
567  void printInfo() const;
568 
569  void setRelativeRho(int start, int end, EGS_Float rho);
570  void setRelativeRho(EGS_Input *);
571  EGS_Float getRelativeRho(int ireg) const {
572  if (ireg < 0 || ireg >= nreg) {
573  return 1;
574  }
575  if (ireg < nbase) {
576  return g->getRelativeRho(ireg);
577  }
578  int jg = (ireg - nbase)/nmax;
579  return geometries[jg]->getRelativeRho(ireg - nbase - jg*nmax);
580  };
581 
582  void setBScaling(int start, int end, EGS_Float bf);
583  void setBScaling(EGS_Input *);
584  EGS_Float getBScaling(int ireg) const {
585  if (ireg < 0 || ireg >= nreg) {
586  return 1;
587  }
588  if (ireg < nbase) {
589  return g->getBScaling(ireg);
590  }
591  int jg = (ireg - nbase)/nmax;
592  return geometries[jg]->getBScaling(ireg - nbase - jg*nmax);
593  };
594 
595  virtual void getLabelRegions(const string &str, vector<int> &regs);
596 
597 protected:
598 
601  int n_in;
602  int nbase,
603  nmax;
605  static string type;
606 
609  int *local_start;
610 
617  void setMedia(EGS_Input *,int,const int *);
618 
619 private:
620 
621  void setPropertyError(const char *funcname) {
622  egsFatal("EGS_EnvelopeGeometry::%s: don't use this method\n Define "
623  "properties in the constituent geometries instead\n",
624  funcname);
625  };
626 
627 };
628 
629 
630 struct EnvelopeAux;
631 
640 class EGS_ENVELOPEG_EXPORT EGS_FastEnvelope : public EGS_BaseGeometry {
641 
642 public:
643 
645  const vector<EnvelopeAux *> &fgeoms, const string &Name = "",
646  int newindexing=false);
647 
648  ~EGS_FastEnvelope();
649 
650  bool isRealRegion(int ireg) const {
651  if (ireg < 0 || ireg >= nreg) {
652  return false;
653  }
654  if (ireg < nbase) {
655  return g->isRealRegion(ireg);
656  }
657  int jg, ilocal;
658  if (new_indexing) {
659  jg = reg_to_inscr[ireg-nbase];
660  ilocal = ireg - local_start[jg];
661  }
662  else {
663  jg = (ireg - nbase)/nmax;
664  ilocal = ireg - nbase - jg*nmax;
665  }
666  return geometries[jg]->isRealRegion(ilocal);
667  };
668 
669  bool isInside(const EGS_Vector &x) {
670  return g->isInside(x);
671  };
672 
673  int isWhere(const EGS_Vector &x) {
674  int ireg = g->isWhere(x);
675  if (ireg < 0 || n_start[ireg] < 0) {
676  return ireg;
677  }
678  for (int jj=n_start[ireg]; jj<n_start[ireg+1]; jj++) {
679  int j = glist[jj];
680  int i = geometries[j]->isWhere(x);
681  if (i >= 0) {
682  return nbase + nmax*j + i;
683  }
684  }
685  return ireg;
686  };
687 
688  int inside(const EGS_Vector &x) {
689  return isWhere(x);
690  };
691 
692  int medium(int ireg) const {
693  if (ireg < nbase) {
694  return g->medium(ireg);
695  }
696  int jg = (ireg - nbase)/nmax;
697  int ilocal = ireg - nbase - jg*nmax;
698  return geometries[jg]->medium(ilocal);
699  };
700 
701  int computeIntersections(int ireg, int n, const EGS_Vector &X,
702  const EGS_Vector &u, EGS_GeometryIntersections *isections) {
703  if (n < 1) {
704  return -1;
705  }
706  int ifirst = 0;
707  EGS_Float t, ttot = 0;
708  EGS_Vector x(X);
709  int imed;
710  //egsInformation("computeIntersections: ireg=%d x=(%g,%g,%g) "
711  // "u=(%g,%g,%g)\n",ireg,x.x,x.y,x.z,u.x,u.y,u.z);
712  if (ireg < 0) {
713  t = veryFar;
714  ireg = howfar(ireg,x,u,t,&imed);
715  if (ireg < 0) {
716  return 0;
717  }
718  isections[0].t = t;
719  isections[0].rhof = 1;
720  isections[0].ireg = -1;
721  isections[0].imed = -1;
722  ttot = t;
723  ++ifirst;
724  x += u*t;
725  //egsInformation("entered after t=%g in ireg=%d at x=(%g,%g,%g)\n",
726  // t,ireg,x.x,x.y,x.z);
727  }
728  else {
729  imed = medium(ireg);
730  }
731 
732 
733  int j = ifirst;
734  int ij = -1, ig=0;
735  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
736  if (loopCount == loopMax) {
737  egsFatal("EGS_FastEnvelope::computeIntersections: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
738  return -1;
739  }
740  //egsInformation("in loop: j=%d ireg=%d imed=%d x=(%g,%g,%g)\n",
741  // j,ireg,imed,x.x,x.y,x.z);
742  if (ireg >= nbase) {
743  ig = (ireg - nbase)/nmax;
744  ij = ireg - nbase - ig*nmax;
745  }
746  isections[j].imed = imed;
747  isections[j].ireg = ireg;
748  isections[j].rhof = getRelativeRho(ireg);
749  if (ireg < nbase) { // in one of the regions of the base geometry
750  t = veryFar;
751  int ibase = g->howfar(ireg,x,u,t,&imed);
752  //egsInformation("In base geometry: t=%g inew=%d\n",t,ibase);
753  ij = -1;
754  for (int ii=n_start[ireg]; ii<n_start[ireg+1]; ii++) {
755  int i = glist[ii];
756  int ireg_i = geometries[i]->howfar(-1,x,u,t,&imed);
757  if (ireg_i >= 0) {
758  ij = ireg_i;
759  ig = i;
760  }
761  }
762  ttot += t;
763  isections[j++].t = ttot;
764  //egsInformation("after inscribed loop: t=%g ij=%d ig=%d"
765  // " ttot=%g\n",t,ij,ig,ttot);
766  if (ij < 0) {
767  ireg = ibase;
768  }
769  else {
770  ireg = nbase + ig*nmax + ij;
771  }
772  if (ireg < 0) {
773  return j;
774  }
775  if (j >= n) {
776  return -1;
777  }
778  x += u*t;
779  }
780  else {
781  int iadd = nbase + ig*nmax;
782  int nsec = geometries[ig]->computeIntersections(ij,n-j,
783  x,u,&isections[j]);
784  //egsInformation("In inscribed %d: got %d intersections\n",ig,
785  // nsec);
786  int nm = nsec >= 0 ? nsec+j : n;
787  for (int i=j; i<nm; i++) {
788  isections[i].ireg += iadd;
789  isections[i].t += ttot;
790  }
791  //egsInformation("last intersection: %g\n",isections[nm-1].t);
792  if (nsec < 0) {
793  return nsec;
794  }
795  j += nsec;
796  if (j >= n) {
797  return -1;
798  }
799  t = isections[j-1].t - ttot;
800  x += u*t;
801  ttot = isections[j-1].t;
802  ireg = g->isWhere(x);
803  //egsInformation("new region: %d\n",ireg);
804  if (ireg < 0) {
805  return j;
806  }
807  imed = g->medium(ireg);
808  }
809  }
810  return -1;
811  }
812 
813  EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
814  const EGS_Vector &u) {
815  if (ireg < 0) {
816  return 0;
817  }
818  EGS_Float d;
819  if (ireg < nbase) {
820  d = g->howfarToOutside(ireg,x,u);
821  }
822  else if (g->regions() == 1) {
823  d = g->howfarToOutside(0,x,u);
824  }
825  else {
826  int ir = g->isWhere(x);
827  d = g->howfarToOutside(ir,x,u);
828  }
829  return d;
830  };
831 
832  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
833  EGS_Float &t, int *newmed = 0, EGS_Vector *normal = 0) {
834  if (ireg >= 0) {
835  // inside.
836  if (ireg < nbase) {
837  // in one of the regions of the base geometry
838  // check if we hit a boundary in the base geometry.
839  // if we do, newmed and normal get set accordingly.
840  int ibase = g->howfar(ireg,x,u,t,newmed,normal);
841  int ij = -1, jg;
842  // check if we will enter any of the inscribed geometries
843  // before entering a new region in the base geometry.
844  for (int jj=n_start[ireg]; jj<n_start[ireg+1]; jj++) {
845  int j = glist[jj];
846  int ireg_j =
847  geometries[j]->howfar(-1,x,u,t,newmed,normal);
848  if (ireg_j >= 0) {
849  // we do. remember the inscribed geometry index
850  // and local region
851  ij = ireg_j;
852  jg = j;
853  }
854  }
855  if (ij < 0) {
856  return ibase;
857  }
858  // ij<0 implies that we have not hit any of the
859  // inscribed geometries => return the base geometry index.
860  // ij>=0 implies that we entered inscribed geometry
861  // jg in its local region ij.
862  return nbase + jg*nmax + ij;
863  }
864  // if here, we are in an inscribed geometry.
865  // calculate its index (jg) and its local region (ilocal).
866  int jg = (ireg - nbase)/nmax;
867  int ilocal = ireg - nbase - jg*nmax;
868  // and then check if we will hit a boundary in this geometry.
869  int inew = geometries[jg]->howfar(ilocal,x,u,t,newmed,normal);
870  if (inew >= 0) {
871  return nbase + jg*nmax + inew;
872  }
873  // inew >= 0 implies that we either stay in the same
874  // region (inew=ilocal) or we entered a new region
875  // (inew!=ilocal), which is still inside the inscribed geometry
876  // inew<0 implies that we have exited the inscribed geometry
877  // => check to see in which base geometry region we are.
878  inew = g->isWhere(x+u*t);
879  if (inew >= 0 && newmed) {
880  *newmed = g->medium(inew);
881  }
882  return inew;
883  }
884  // if here, we are outside the base geometry.
885  // check to see if we will enter.
886  int ienter = g->howfar(ireg,x,u,t,newmed,normal);
887  if (ienter >= 0) {
888  // yes, we do. see if we are already inside of one of the
889  // inscribed geometries.
890  //if( n_start[ienter] < 0 ) return ienter;
891  for (int jj=n_start[ienter]; jj<n_start[ienter+1]; jj++) {
892  int j = glist[jj];
893  int i = geometries[j]->isWhere(x+u*t);
894  if (i >= 0) {
895  // yes, we are.
896  if (newmed) {
897  *newmed = geometries[j]->medium(i);
898  }
899  return nbase + nmax*j + i;
900  }
901  }
902  }
903  return ienter;
904  };
905 
906  EGS_Float hownear(int ireg, const EGS_Vector &x) {
907  if (ireg >= 0) {
908  EGS_Float tmin;
909  if (ireg < nbase) { // in one of the regions of the base geom.
910  tmin = g->hownear(ireg,x);
911  if (tmin <= 0) {
912  return tmin;
913  }
914  for (int jj=n_start[ireg]; jj<n_start[ireg+1]; jj++) {
915  int j = glist[jj];
916  EGS_Float tj = geometries[j]->hownear(-1,x);
917  if (tj < tmin) {
918  tmin = tj;
919  if (tmin <= 0) {
920  return tmin;
921  }
922  }
923  }
924  return tmin;
925  }
926  int jg = (ireg - nbase)/nmax;
927  int ilocal = ireg - nbase - jg*nmax;
928  return geometries[jg]->hownear(ilocal,x);
929  }
930  return g->hownear(ireg,x);
931  };
932 
933  bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
934  if (ireg >= 0 && ireg < nreg) {
935  if (ireg < nbase) {
936  return g->hasBooleanProperty(ireg,prop);
937  }
938  int jg = (ireg - nbase)/nmax;
939  int ilocal = ireg - nbase - jg*nmax;
940  return geometries[jg]->hasBooleanProperty(ilocal,prop);
941  }
942  return false;
943  };
944  void setBooleanProperty(EGS_BPType) {
945  setPropertyError("setBooleanProperty()");
946  };
947  void addBooleanProperty(int) {
948  setPropertyError("addBooleanProperty()");
949  };
950  void setBooleanProperty(EGS_BPType,int,int,int step=1) {
951  setPropertyError("setBooleanProperty()");
952  };
953  void addBooleanProperty(int,int,int,int step=1) {
954  setPropertyError("addBooleanProperty()");
955  };
956 
957  int getMaxStep() const {
958  int nstep = g->getMaxStep();
959  for (int j=0; j<n_in; ++j) {
960  nstep += geometries[j]->getMaxStep();
961  }
962  return nstep+n_in;
963  };
964 
965  const string &getType() const {
966  return type;
967  };
968 
969  void printInfo() const;
970 
971  void setRelativeRho(int start, int end, EGS_Float rho);
972  void setRelativeRho(EGS_Input *);
973  EGS_Float getRelativeRho(int ireg) const {
974  if (ireg < 0 || ireg >= nreg) {
975  return 1;
976  }
977  if (ireg < nbase) {
978  return g->getRelativeRho(ireg);
979  }
980  int jg = (ireg - nbase)/nmax;
981  return geometries[jg]->getRelativeRho(ireg - nbase - jg*nmax);
982  };
983 
984  void setBScaling(int start, int end, EGS_Float bf);
985  void setBScaling(EGS_Input *);
986  EGS_Float getBScaling(int ireg) const {
987  if (ireg < 0 || ireg >= nreg) {
988  return 1;
989  }
990  if (ireg < nbase) {
991  return g->getBScaling(ireg);
992  }
993  int jg = (ireg - nbase)/nmax;
994  return geometries[jg]->getBScaling(ireg - nbase - jg*nmax);
995  };
996 
997  virtual void getLabelRegions(const string &str, vector<int> &regs);
998 
999 protected:
1000 
1003  int n_in;
1004  int nbase,
1005  nmax;
1007  int *glist;
1008  int *n_start;
1009  static string type;
1010 
1014 
1021  void setMedia(EGS_Input *,int,const int *);
1022 
1023 private:
1024 
1025  void setPropertyError(const char *funcname) {
1026  egsFatal("EGS_FastEnvelope::%s: don't use this method\n Define "
1027  "properties in the constituent geometries instead\n",
1028  funcname);
1029  };
1030 
1031 
1032 };
1033 
1034 #endif
EGS_BaseGeometry * g
The envelope geometry.
virtual const string & getType() const =0
Get the geometry type.
EGS_BaseGeometry ** geometries
The inscribed geometries.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
int * local_start
First region for each inscribed geometry.
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:95
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const
Is the boolean property prop set for region ireg ?
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
An envelope geometry class.
EGS_Float t
distance to next region boundary
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
EGS_Float rhof
relative mass density in that region
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
Definition: egs_vector.h:56
static string type
Geometry type.
Global egspp functions header file.
EGS_BaseGeometry * g
The envelope geometry.
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.
const EGS_Float veryFar
A very large float.
EGS_I32 imed
medium index
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
int n_in
Number of inscribed geometries.
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
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.
int * reg_to_inscr
Region to inscribed geometry conversion.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
bool new_indexing
If true, use new indexing style.
virtual int computeIntersections(int ireg, int n, const EGS_Vector &x, const EGS_Vector &u, EGS_GeometryIntersections *isections)
Calculates intersection distances to region boundaries.
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
EGS_I32 ireg
region index
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 envelope geometry class.
bool new_indexing
If true, use new indexing style.
EGS_BaseGeometry class header file.
int n_in
Number of inscribed geometries.
int * local_start
First region for each inscribed geometry.
int * reg_to_inscr
Region to inscribed geometry conversion.
EGS_BaseGeometry ** geometries
The inscribed geometries.
virtual void addBooleanProperty(int bit)
Add a boolean property for the entire geometry by setting the bit&#39;th bit.
static string type
Geometry type.
virtual void setBooleanProperty(EGS_BPType prop)
Set the boolean properties of the entire geometry to prop.