45 #include "../egs_gtransformed/egs_gtransformed.h"
47 #ifndef EGS_LATTICE_GEOMETRY_
48 #define EGS_LATTICE_GEOMETRY_
52 #ifdef BUILD_LATTICE_DLL
53 #define EGS_LATTICE_EXPORT __declspec(dllexport)
55 #define EGS_LATTICE_EXPORT __declspec(dllimport)
57 #define EGS_LATTICE_LOCAL
61 #ifdef HAVE_VISIBILITY
62 #define EGS_LATTICE_EXPORT __attribute__ ((visibility ("default")))
63 #define EGS_LATTICE_LOCAL __attribute__ ((visibility ("hidden")))
65 #define EGS_LATTICE_EXPORT
66 #define EGS_LATTICE_LOCAL
158 EGS_Float y, EGS_Float z,
const string &Name =
"");
164 int(round(x.
z/c))*c);
169 return base->computeIntersections(ireg,n,x,u,isections);
174 return base->isRealRegion(ireg);
176 return sub->isRealRegion(ireg - base->regions());
180 return base->isInside(x);
184 int temp = base->isWhere(x);
186 sub->setTransformation(closestPoint(x));
187 if (sub->isInside(x)) {
188 return sub->isWhere(x) + base->regions();
198 int medium(
int ireg)
const {
199 if (ireg >= base->regions()) {
200 return sub->medium(ireg-base->regions());
202 return base->medium(ireg);
211 if (ireg >= base->regions()) {
212 return base->howfarToOutside(ind,x,u);
214 return base->howfarToOutside(ireg,x,u);
218 EGS_Float &t,
int *newmed=0,
EGS_Vector *normal=0) {
225 if (ireg >= base->regions()) {
226 sub->setTransformation(closestPoint(x));
230 int tempReg = sub->howfar(ireg-base->regions(),x,u,tempT,newmed,normal);
237 newX = x + (newX * tempT);
239 if (base->isWhere(newX) != ind) {
242 tempReg = base->howfar(ind,x,u,t,0,normal);
243 if (newmed && tempReg >= 0) {
244 *newmed = base->medium(tempReg);
256 *newmed = base->medium(ind);
268 return tempReg+base->regions();
270 else if (ireg == ind) {
273 base->howfar(ireg,x,u,tempT);
274 EGS_Float max = tempT;
277 EGS_Float min = max, minX = max, minY = max, minZ = max;
280 EGS_Vector xInt = unit*(a/4.0/fabs(unit.
x)), yInt = unit*(b/4.0/fabs(unit.
y)), zInt = unit*(c/4.0/fabs(unit.
z));
285 EGS_Float max2 = max*max;
288 while ((x-x0).length2() < max2) {
290 tempP = closestPoint(x0);
291 sub->setTransformation(tempP);
292 if (sub->howfar(-1,x,u,tempT)+1)
293 if (tempT < max && tempT > 0) {
302 while ((x-x0).length2() < max2) {
304 tempP = closestPoint(x0);
305 sub->setTransformation(tempP);
306 if (sub->howfar(-1,x,u,tempT)+1)
307 if (tempT < max && tempT > 0) {
316 while ((x-x0).length2() < max2) {
318 tempP = closestPoint(x0);
319 sub->setTransformation(tempP);
320 if (sub->howfar(-1,x,u,tempT)+1)
321 if (tempT < max && tempT > 0) {
342 sub->setTransformation(closestPoint(x+unit*tempT));
343 sub->howfar(-1,x,u,tempT);
344 if (tempT < min && tempT > 0) {
346 finalP = closestPoint(x+unit*tempT);
352 sub->setTransformation(finalP);
353 int tempReg = sub->howfar(-1,x,u,tempT,newmed,normal);
355 if (newmed && tempReg >= 0) {
356 *newmed = sub->medium(tempReg);
364 return tempReg+base->regions();
369 int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
374 if (newmed && tempReg >= 0) {
375 *newmed = base->medium(tempReg);
383 int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
386 if (tempReg == ind) {
390 newX = x + (newX * t);
392 sub->setTransformation(closestPoint(newX));
393 int newReg = sub->isWhere(newX);
395 if (newmed && newReg >= 0) {
396 *newmed = sub->medium(newReg);
402 return newReg + base->regions();
406 if (newmed && tempReg >= 0) {
407 *newmed = base->medium(tempReg);
418 EGS_Float temp, dist = base->hownear(ind,x);
419 if (ireg >= base->regions()) {
420 sub->setTransformation(closestPoint(x));
421 temp = sub->hownear(ireg-base->regions(),x);
422 dist=(temp<dist)?temp:dist;
424 else if (ireg == ind) {
426 EGS_Float xa = floor(x.
x/a)*a, yb = floor(x.
y/b)*b, zc = floor(x.
z/c)*c;
436 for (
int i = 0; i < 8; i++) {
437 sub->setTransformation(x0[i]);
438 temp = sub->hownear(-1,x);
445 dist = base->hownear(ireg,x);
451 return base->regions() + sub->regions();
458 const string &
getType()
const {
467 return base->hasBooleanProperty(ireg, prop);
469 return sub->hasBooleanProperty(ireg - base->regions(), prop);
472 base->setBooleanProperty(prop);
475 base->addBooleanProperty(bit);
478 base->setBooleanProperty(prop,start,end,step);
481 base->addBooleanProperty(bit,start,end,step);
486 return base->getRelativeRho(ireg);
488 return sub->getRelativeRho(ireg - base->regions());
493 void setBScaling(
int start,
int end, EGS_Float rho);
497 return base->getBScaling(ireg);
499 return sub->getBScaling(ireg - base->regions());
528 EGS_Float i1,j1,k1,i2,j2,j3,k3,j4;
529 EGS_Float r3a = 2.0*gap;
531 i1 = a*(round(x.
x/a));
532 j1 = r3a*(round(x.
y/r3a));
533 k1 = r3a*(round(x.
z/r3a));
535 i2 = a*(0.5+round(x.
x/a-0.5));
536 j2 = r3a*(0.5+round(x.
y/r3a-0.5));
538 j3 = r3a*(0.25+round(x.
y/r3a-0.25));
539 k3 = r3a*(0.5+round(x.
z/r3a-0.5));
541 j4 = r3a*(-0.25+round(x.
y/r3a+0.25));
548 d[0] = (x-p1).length2();
549 d[1] = (x-p2).length2();
550 d[2] = (x-p3).length2();
551 d[3] = (x-p4).length2();
553 if (d[0] < d[1] && d[0] < d[2] && d[0] < d[3]) {
556 if (d[1] < d[2] && d[1] < d[3]) {
567 return base->computeIntersections(ireg,n,x,u,isections);
572 return base->isRealRegion(ireg);
574 return sub->isRealRegion(ireg - base->regions());
578 return base->isInside(x);
582 int temp = base->isWhere(x);
584 sub->setTransformation(closestPoint(x));
585 if (sub->isInside(x)) {
586 return sub->isWhere(x) + base->regions();
596 int medium(
int ireg)
const {
597 if (ireg >= base->regions()) {
598 return sub->medium(ireg-base->regions());
600 return base->medium(ireg);
609 if (ireg >= base->regions()) {
610 return base->howfarToOutside(ind,x,u);
612 return base->howfarToOutside(ireg,x,u);
616 EGS_Float &t,
int *newmed=0,
EGS_Vector *normal=0) {
623 if (ireg >= base->regions()) {
624 sub->setTransformation(closestPoint(x));
628 int tempReg = sub->howfar(ireg-base->regions(),x,u,tempT,newmed,normal);
635 newX = x + (newX * tempT);
637 if (base->isWhere(newX) != ind) {
640 tempReg = base->howfar(ind,x,u,t,0,normal);
641 if (newmed && tempReg >= 0) {
642 *newmed = base->medium(tempReg);
654 *newmed = base->medium(ind);
666 return tempReg+base->regions();
668 else if (ireg == ind) {
671 base->howfar(ireg,x,u,tempT);
672 EGS_Float max = tempT;
675 EGS_Float min = max, minX = max, minY = max, minZ = max;
678 EGS_Vector xInt = unit*(a/8.0/fabs(unit.
x)), yInt = unit*(gap/8.0/fabs(unit.
y)), zInt = unit*(gap/8.0/fabs(unit.
z));
683 EGS_Float max2 = max*max;
686 while ((x-x0).length2() < max2) {
688 tempP = closestPoint(x0);
689 sub->setTransformation(tempP);
690 if (sub->howfar(-1,x,u,tempT)+1)
691 if (tempT < max && tempT > 0) {
700 while ((x-x0).length2() < max2) {
702 tempP = closestPoint(x0);
703 sub->setTransformation(tempP);
704 if (sub->howfar(-1,x,u,tempT)+1)
705 if (tempT < max && tempT > 0) {
714 while ((x-x0).length2() < max2) {
716 tempP = closestPoint(x0);
717 sub->setTransformation(tempP);
718 if (sub->howfar(-1,x,u,tempT)+1)
719 if (tempT < max && tempT > 0) {
740 sub->setTransformation(closestPoint(x+unit*tempT));
741 sub->howfar(-1,x,u,tempT);
742 if (tempT < min && tempT > 0) {
744 finalP = closestPoint(x+unit*tempT);
750 sub->setTransformation(finalP);
751 int tempReg = sub->howfar(-1,x,u,tempT,newmed,normal);
753 if (newmed && tempReg >= 0) {
754 *newmed = sub->medium(tempReg);
762 return tempReg+base->regions();
767 int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
772 if (newmed && tempReg >= 0) {
773 *newmed = base->medium(tempReg);
781 int tempReg = base->howfar(ireg,x,u,t,newmed,normal);
784 if (tempReg == ind) {
788 newX = x + (newX * t);
790 sub->setTransformation(closestPoint(newX));
791 int newReg = sub->isWhere(newX);
793 if (newmed && newReg >= 0) {
794 *newmed = sub->medium(newReg);
800 return newReg + base->regions();
804 if (newmed && tempReg >= 0) {
805 *newmed = base->medium(tempReg);
816 EGS_Float temp, dist = base->hownear(ireg>=base->regions()?ind:ireg,x);
817 if (ireg >= base->regions()) {
818 sub->setTransformation(closestPoint(x));
819 temp = sub->hownear(ireg-base->regions(),x);
820 dist=(temp<dist)?temp:dist;
822 else if (ireg == ind) {
825 EGS_Float r3a = 2.0*gap, hgap = gap/2.0, ha = a/2.0;
829 i = a*(floor(x.
x/a));
830 j = r3a*(round(x.
y/r3a));
831 k = r3a*(round(x.
z/r3a));
835 if (x.
y<=j && x.
z<=k)
845 else if (x.
y>j && x.
z<=k)
855 else if (x.
y<j && x.
z>k)
875 for (
int i = 0; i < 8; i++) {
876 sub->setTransformation(x0[i]);
877 temp = sub->hownear(-1,x);
887 return base->regions() + sub->regions();
894 const string &
getType()
const {
903 return base->hasBooleanProperty(ireg, prop);
905 return sub->hasBooleanProperty(ireg - base->regions(), prop);
908 base->setBooleanProperty(prop);
911 base->addBooleanProperty(bit);
914 base->setBooleanProperty(prop,start,end,step);
917 base->addBooleanProperty(bit,start,end,step);
922 return base->getRelativeRho(ireg);
924 return sub->getRelativeRho(ireg - base->regions());
929 void setBScaling(
int start,
int end, EGS_Float rho);
933 return base->getBScaling(ireg);
935 return sub->getBScaling(ireg - base->regions());
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const
Is the boolean property prop set for region ireg ?
int regions() const
Returns the number of local regions in this geometry.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
vector< EGS_Float > d
Don't redefine 4 dists in closestPoint each invocation.
EGS_Float c
The center-to-center distance along x, y, and z.
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
EGS_TransformedGeometry * sub
The sub geometry that could appear within base.
A class representing 3D vectors.
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.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
A Bravais, cubic, and hexagonal lattice geometryA geometry which embeds a lattice of one geometry (na...
EGS_BaseGeometry * base
The geometry within which the sub geometry appears.
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.
virtual void getLabelRegions(const string &str, vector< int > ®s)
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.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
virtual int computeIntersections(int ireg, int n, const EGS_Vector &x, const EGS_Vector &u, EGS_GeometryIntersections *isections)
Calculates intersection distances to region boundaries.
EGS_Float a
The center-to-center distance to the nearest 12 neighbours.
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.
int ind
The region in base geom where we could encounter sub geom.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
EGS_Float gap
Translating the hexagonal lattice to xyz coordinate,.
int ind
The region in base geom where we could encounter sub geom.
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
string type
The geometry type.
EGS_TransformedGeometry * sub
The sub geometry that could appear within base.
string type
The geometry type.
EGS_BaseGeometry class header file.
EGS_BaseGeometry * base
The geometry within which the sub geometry appears.
int maxStep
The maximum number of steps.
virtual void addBooleanProperty(int bit)
Add a boolean property for the entire geometry by setting the bit'th bit.
virtual void setBooleanProperty(EGS_BPType prop)
Set the boolean properties of the entire geometry to prop.
int maxStep
The maximum number of steps.