Line data    Source code

       1             : #include "crpropa/module/Boundary.h"
       2             : #include "crpropa/Units.h"
       3             : 
       4             : #include <sstream>
       5             : 
       6             : namespace crpropa {
       7             : 
       8           0 : PeriodicBox::PeriodicBox() :
       9           0 :                 origin(Vector3d(0, 0, 0)), size(Vector3d(0, 0, 0)) {
      10           0 : }
      11             : 
      12           2 : PeriodicBox::PeriodicBox(Vector3d o, Vector3d s) :
      13           2 :                 origin(o), size(s) {
      14           2 : }
      15             : 
      16           2 : void PeriodicBox::process(Candidate *c) const {
      17           2 :         Vector3d pos = c->current.getPosition();
      18           2 :         Vector3d n = ((pos - origin) / size).floor();
      19             : 
      20           2 :         if ((n.x == 0) and (n.y == 0) and (n.z == 0))
      21             :                 return; // do nothing if candidate is inside the box
      22             : 
      23           2 :         c->current.setPosition(pos - n * size);
      24           2 :         c->previous.setPosition(c->previous.getPosition() - n * size);
      25           2 :         c->source.setPosition(c->source.getPosition() - n * size);
      26           2 :         c->created.setPosition(c->created.getPosition() - n * size);
      27             : }
      28             : 
      29           0 : void PeriodicBox::setOrigin(Vector3d o) {
      30             :         origin = o;
      31           0 : }
      32           0 : void PeriodicBox::setSize(Vector3d s) {
      33             :         size = s;
      34           0 : }
      35             : 
      36           0 : std::string PeriodicBox::getDescription() const {
      37           0 :         std::stringstream s;
      38           0 :         s << "Periodic box: origin " << origin / Mpc << " Mpc, ";
      39           0 :         s << "size " << size / Mpc << " Mpc";
      40           0 :         return s.str();
      41           0 : }
      42             : 
      43           0 : ReflectiveBox::ReflectiveBox() :
      44           0 :                 origin(Vector3d(0, 0, 0)), size(Vector3d(0, 0, 0)) {
      45           0 : }
      46             : 
      47           1 : ReflectiveBox::ReflectiveBox(Vector3d o, Vector3d s) :
      48           1 :                 origin(o), size(s) {
      49           1 : }
      50             : 
      51           1 : void ReflectiveBox::process(Candidate *c) const {
      52           1 :         Vector3d cur = (c->current.getPosition() - origin) / size; // current position in cell units
      53           1 :         Vector3d n = cur.floor();
      54             : 
      55           1 :         if ((n.x == 0) and (n.y == 0) and (n.z == 0))
      56             :                 return; // do nothing if candidate is inside the box
      57             : 
      58             :         // flip direction
      59           1 :         Vector3d nReflect(pow(-1, n.x), pow(-1, n.y), pow(-1, n.z));
      60           1 :         c->current.setDirection(c->current.getDirection() * nReflect);
      61           1 :         c->previous.setDirection(c->previous.getDirection() * nReflect);
      62           1 :         c->created.setDirection(c->created.getDirection() * nReflect);
      63           1 :         c->source.setDirection(c->source.getDirection() * nReflect);
      64             : 
      65           1 :         Vector3d src = (c->source.getPosition() - origin) / size; // initial position in cell units
      66           1 :         Vector3d cre = (c->created.getPosition() - origin) / size; // initial position in cell units
      67           1 :         Vector3d prv = (c->previous.getPosition() - origin) / size; // previous position in cell units
      68             : 
      69             :         // repeatedly translate until the current position is inside the cell
      70           1 :         while ((cur.x < 0) or (cur.x > 1)) {
      71           0 :                 double t = 2 * (cur.x > 1);
      72           0 :                 src.x = t - src.x;
      73           0 :                 cre.x = t - cre.x;
      74           0 :                 prv.x = t - prv.x;
      75           0 :                 cur.x = t - cur.x;
      76             :         }
      77           1 :         while ((cur.y < 0) or (cur.y > 1)) {
      78           0 :                 double t = 2 * (cur.y > 1);
      79           0 :                 src.y = t - src.y;
      80           0 :                 cre.y = t - cre.y;
      81           0 :                 prv.y = t - prv.y;
      82           0 :                 cur.y = t - cur.y;
      83             :         }
      84           2 :         while ((cur.z < 0) or (cur.z > 1)) {
      85           1 :                 double t = 2 * (cur.z > 1);
      86           1 :                 src.z = t - src.z;
      87           1 :                 cre.z = t - cre.z;
      88           1 :                 prv.z = t - prv.z;
      89           1 :                 cur.z = t - cur.z;
      90             :         }
      91             : 
      92           1 :         c->current.setPosition(cur * size + origin);
      93           1 :         c->source.setPosition(src * size + origin);
      94           1 :         c->created.setPosition(cre * size + origin);
      95           1 :         c->previous.setPosition(prv * size + origin);
      96             : }
      97             : 
      98           0 : void ReflectiveBox::setOrigin(Vector3d o) {
      99             :         origin = o;
     100           0 : }
     101           0 : void ReflectiveBox::setSize(Vector3d s) {
     102             :         size = s;
     103           0 : }
     104             : 
     105           0 : std::string ReflectiveBox::getDescription() const {
     106           0 :         std::stringstream s;
     107           0 :         s << "Reflective box: origin " << origin / Mpc << " Mpc, ";
     108           0 :         s << "size " << size / Mpc << " Mpc";
     109           0 :         return s.str();
     110           0 : }
     111             : 
     112           0 : CubicBoundary::CubicBoundary() :
     113           0 :                 origin(Vector3d(0, 0, 0)), size(0), limitStep(true), margin(0.1 * kpc) {
     114           0 : }
     115             : 
     116           4 : CubicBoundary::CubicBoundary(Vector3d o, double s) :
     117           4 :                 origin(o), size(s), limitStep(true), margin(0.1 * kpc) {
     118           4 : }
     119             : 
     120           4 : void CubicBoundary::process(Candidate *c) const {
     121           4 :         Vector3d r = c->current.getPosition() - origin;
     122             :         double lo = r.min();
     123             :         double hi = r.max();
     124           4 :         if ((lo <= 0) or (hi >= size)) {
     125           1 :                 reject(c);
     126             :         }
     127           4 :         if (limitStep) {
     128           4 :                 c->limitNextStep(lo + margin);
     129           4 :                 c->limitNextStep(size - hi + margin);
     130             :         }
     131           4 : }
     132             : 
     133           0 : void CubicBoundary::setOrigin(Vector3d o) {
     134             :         origin = o;
     135           0 : }
     136           0 : void CubicBoundary::setSize(double s) {
     137           0 :         size = s;
     138           0 : }
     139           2 : void CubicBoundary::setMargin(double m) {
     140           2 :         margin = m;
     141           2 : }
     142           2 : void CubicBoundary::setLimitStep(bool b) {
     143           2 :         limitStep = b;
     144           2 : }
     145             : 
     146           0 : std::string CubicBoundary::getDescription() const {
     147           0 :         std::stringstream s;
     148           0 :         s << "Cubic Boundary: origin " << origin / Mpc << " Mpc, ";
     149           0 :         s << "size " << size / Mpc << " Mpc, ";
     150           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     151           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     152           0 :         if (rejectAction.valid())
     153           0 :                 s << ", Action: " << rejectAction->getDescription();
     154           0 :         return s.str();
     155           0 : }
     156             : 
     157           0 : SphericalBoundary::SphericalBoundary() :
     158           0 :                 center(Vector3d(0, 0, 0)), radius(0), limitStep(true), margin(0.1 * kpc) {
     159           0 : }
     160             : 
     161           3 : SphericalBoundary::SphericalBoundary(Vector3d c, double r) :
     162           3 :                 center(c), radius(r), limitStep(true), margin(0.1 * kpc) {
     163           3 : }
     164             : 
     165           3 : void SphericalBoundary::process(Candidate *c) const {
     166           3 :         double d = (c->current.getPosition() - center).getR();
     167           3 :         if (d >= radius) {
     168           1 :                 reject(c);
     169             :         }
     170           3 :         if (limitStep)
     171           3 :                 c->limitNextStep(radius - d + margin);
     172           3 : }
     173             : 
     174           0 : void SphericalBoundary::setCenter(Vector3d c) {
     175             :         center = c;
     176           0 : }
     177           0 : void SphericalBoundary::setRadius(double r) {
     178           0 :         radius = r;
     179           0 : }
     180           1 : void SphericalBoundary::setMargin(double m) {
     181           1 :         margin = m;
     182           1 : }
     183           1 : void SphericalBoundary::setLimitStep(bool b) {
     184           1 :         limitStep = b;
     185           1 : }
     186             : 
     187           0 : std::string SphericalBoundary::getDescription() const {
     188           0 :         std::stringstream s;
     189           0 :         s << "Spherical Boundary: radius " << radius / Mpc << " Mpc, ";
     190           0 :         s << "around " << center / Mpc << " Mpc, ";
     191           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     192           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     193           0 :         if (rejectAction.valid())
     194           0 :                 s << ", Action: " << rejectAction->getDescription();
     195           0 :         return s.str();
     196           0 : }
     197             : 
     198           0 : EllipsoidalBoundary::EllipsoidalBoundary() :
     199             :                 focalPoint1(Vector3d(0, 0, 0)), focalPoint2(Vector3d(0, 0, 0)),
     200           0 :                 majorAxis(0), limitStep(true), margin(0.1 * kpc) {
     201           0 : }
     202             : 
     203           3 : EllipsoidalBoundary::EllipsoidalBoundary(Vector3d f1, Vector3d f2, double a) :
     204           3 :                 focalPoint1(f1), focalPoint2(f2), majorAxis(a), limitStep(true),
     205           3 :                 margin(0.1 * kpc) {
     206           3 : }
     207             : 
     208           3 : void EllipsoidalBoundary::process(Candidate *c) const {
     209           3 :         Vector3d pos = c->current.getPosition();
     210           3 :         double d = pos.getDistanceTo(focalPoint1) + pos.getDistanceTo(focalPoint2);
     211           3 :         if (d >= majorAxis) {
     212           1 :                 reject(c);
     213             :         }
     214           3 :         if (limitStep)
     215           3 :                 c->limitNextStep(majorAxis - d + margin);
     216           3 : }
     217             : 
     218           0 : void EllipsoidalBoundary::setFocalPoints(Vector3d f1, Vector3d f2) {
     219             :         focalPoint1 = f1;
     220             :         focalPoint2 = f2;
     221           0 : }
     222           0 : void EllipsoidalBoundary::setMajorAxis(double a) {
     223           0 :         majorAxis = a;
     224           0 : }
     225           1 : void EllipsoidalBoundary::setMargin(double m) {
     226           1 :         margin = m;
     227           1 : }
     228           1 : void EllipsoidalBoundary::setLimitStep(bool b) {
     229           1 :         limitStep = b;
     230           1 : }
     231             : 
     232           0 : std::string EllipsoidalBoundary::getDescription() const {
     233           0 :         std::stringstream s;
     234           0 :         s << "Ellipsoidal Boundary: F1 = " << focalPoint1 / Mpc << " Mpc, ";
     235           0 :         s << "F2 = " << focalPoint2 / Mpc << " Mpc, ";
     236           0 :         s << "major axis " << majorAxis / Mpc << " Mpc; ";
     237           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     238           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     239           0 :         if (rejectAction.valid())
     240           0 :                 s << ", Action: " << rejectAction->getDescription();
     241           0 :         return s.str();
     242           0 : }
     243             : 
     244           0 : CylindricalBoundary::CylindricalBoundary() :
     245           0 :   origin(Vector3d(0,0,0)), height(0), radius(0), limitStep(false), margin(0) {
     246           0 : }
     247             : 
     248           3 : CylindricalBoundary::CylindricalBoundary(Vector3d o, double h, double r) :
     249           3 :   origin(o), height(h), radius(r), limitStep(false) , margin(0){
     250           3 : }
     251             : 
     252           3 : void CylindricalBoundary::process(Candidate *c) const {
     253           3 :         Vector3d d = c->current.getPosition() - origin;
     254           3 :         double R2 = pow(d.x, 2.)+pow(d.y, 2.);
     255           3 :         double Z = fabs(d.z);
     256           3 :         if ( R2 < pow(radius, 2.) and Z < height/2.) {
     257           2 :                 if(limitStep) {
     258           2 :                         c->limitNextStep(std::min(radius - pow(R2, 0.5), height/2. - Z) + margin);   
     259             :                 }
     260             :           return;
     261             :         }
     262           1 :         reject(c);
     263             : }
     264             : 
     265           0 : void CylindricalBoundary::setOrigin(Vector3d o) {
     266             :         origin = o;
     267           0 : }
     268           0 : void CylindricalBoundary::setHeight(double h) {
     269           0 :         height = h;
     270           0 : }
     271           0 : void CylindricalBoundary::setRadius(double r) {
     272           0 :         radius = r;
     273           0 : }
     274           1 : void CylindricalBoundary::setLimitStep(bool b) {
     275           1 :         limitStep = b;
     276           1 : }
     277           1 : void CylindricalBoundary::setMargin(double m) {
     278           1 :         margin = m;
     279           1 : }
     280             : 
     281             : 
     282           0 : std::string CylindricalBoundary::getDescription() const {
     283           0 :         std::stringstream s;
     284           0 :         s << "Cylindrical Boundary: origin = " << origin / kpc << " kpc, ";
     285           0 :         s << "radius = " << radius / kpc << " kpc, ";
     286           0 :         s << "height" << height / kpc << " kpc; ";
     287           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     288           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     289           0 :         if (rejectAction.valid())
     290           0 :                 s << ", Action: " << rejectAction->getDescription();
     291           0 :         return s.str();
     292           0 : }
     293             : 
     294             : } // namespace crpropa