Line data    Source code

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