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