Line data    Source code

       1              : #include <limits>
       2              : #include <cmath>
       3              : #include "kiss/logger.h"
       4              : #include "crpropa/Geometry.h"
       5              : 
       6              : #include <iostream>
       7              : 
       8              : namespace crpropa {
       9              : 
      10              : // Plane ------------------------------------------------------------------
      11            1 : Plane::Plane(const Vector3d& _x0, const Vector3d& _n) : x0(_x0), n(_n) {
      12            1 : };
      13              : 
      14            0 : Plane::Plane(const Vector3d& _x0, const Vector3d& v1,const Vector3d& v2) : x0(_x0), n(0,0,0) {
      15              :         n = v1.cross(v2);
      16              :         n /= n.getR();
      17            0 : };
      18              : 
      19            2 : double Plane::distance(const Vector3d &x) const {
      20              :         Vector3d dX = x - x0;
      21            2 :         return n.dot(dX);
      22              : };
      23              : 
      24            0 : std::string Plane::getDescription() const {
      25            0 :         std::stringstream ss;
      26              :         ss << "Plane: " << std::endl
      27            0 :            << "   x0: " << x0 << std::endl
      28            0 :            << "    n: " << n << std::endl;
      29            0 :         return ss.str();
      30            0 : };
      31              : 
      32            0 : Vector3d Plane::normal(const Vector3d& point) const {
      33            0 :         return n;
      34              : }
      35              : 
      36              : 
      37              : // Sphere ------------------------------------------------------------------
      38            4 : Sphere::Sphere(const Vector3d& _center, double _radius) : center(_center), radius(_radius) {
      39            4 : };
      40              : 
      41           12 : double Sphere::distance(const Vector3d &point) const {
      42              :         Vector3d dR = point - center;
      43           12 :         return dR.getR() - radius;
      44              : }
      45              : 
      46            0 : Vector3d Sphere::normal(const Vector3d& point) const {
      47              :         Vector3d d = point - center;
      48            0 :         return d.getUnitVector();
      49              : }
      50              : 
      51            0 : std::string Sphere::getDescription() const {
      52            0 :         std::stringstream ss;
      53              :         ss << "Sphere: " << std::endl
      54            0 :                         << "   Center: " << center << std::endl
      55            0 :                         << "   Radius: " << radius << std::endl;
      56            0 :         return ss.str();
      57            0 : };
      58              : 
      59              : 
      60              : // ParaxialBox -------------------------------------------------------------
      61            1 : ParaxialBox::ParaxialBox(const Vector3d& _corner, const Vector3d& _size) : corner(_corner), size(_size) {
      62            1 : };
      63              : 
      64            6 : double ParaxialBox::distance(const Vector3d &point) const {
      65              :         Vector3d X = point - corner - size/2.;
      66              : 
      67              :         // inside the cube
      68            6 :         if ((fabs(X.x) <= size.x/2.) and (fabs(X.y) <= size.y/2.) and (fabs(X.z) <= size.z/2.)) { 
      69              :                 Vector3d Xp = size/2. - X.abs();
      70            3 :                 double d = std::min(Xp.x, std::min(Xp.y, Xp.z));
      71              : 
      72            3 :                 return -1. * d;
      73              :         }
      74              : 
      75            3 :         double a = std::max(0., fabs(X.x) - size.x/2.);
      76            3 :         double b = std::max(0., fabs(X.y) - size.y/2.);
      77            3 :         double c = std::max(0., fabs(X.z) - size.z/2.);
      78              : 
      79            3 :         return sqrt(a*a + b*b +c*c);
      80              : }
      81              : 
      82            0 : Vector3d ParaxialBox::normal(const Vector3d& point) const {
      83              :         Vector3d d = (point - corner).abs();
      84              :         Vector3d d2 = d + size;
      85              :         Vector3d n;
      86              :         double dmin = std::numeric_limits<double>::infinity();
      87            0 :         if (d.x < dmin) {
      88              :                 dmin = d.x;
      89              :                 n = Vector3d(-1, 0, 0);
      90              :         }
      91            0 :         if (d.y < dmin) {
      92              :                 dmin = d.y;
      93              :                 n = Vector3d(0, -1, 0);
      94              :         }
      95            0 :         if (d.z < dmin) {
      96              :                 dmin = d.z;
      97              :                 n = Vector3d(0, 0, -1);
      98              :         }
      99            0 :         if (d2.x < dmin) {
     100              :                 dmin = d2.x;
     101              :                 n = Vector3d(1, 0, 0);
     102              :         }
     103            0 :         if (d2.y < dmin) {
     104              :                 dmin = d2.y;
     105              :                 n = Vector3d(0, 1, 0);
     106              :         }
     107            0 :         if (d2.z < dmin) {
     108              :                 // dmin = d2.z;
     109              :                 n = Vector3d(0, 0, 1);
     110              :         }
     111              : 
     112            0 :         return n;
     113              : }
     114              : 
     115            0 : std::string ParaxialBox::getDescription() const {
     116            0 :         std::stringstream ss;
     117              :         ss << "ParaxialBox: " << std::endl
     118            0 :                  << "   corner: " << corner << std::endl
     119            0 :                  << "     size: " << size << std::endl;
     120            0 :         return ss.str();
     121            0 : };
     122              : 
     123              : 
     124              : } // namespace