Line data    Source code

       1             : #include "crpropa/magneticField/MagneticField.h"
       2             : 
       3             : namespace crpropa {
       4             : 
       5           0 : PeriodicMagneticField::PeriodicMagneticField(ref_ptr<MagneticField> field,
       6           0 :                 const Vector3d &extends) :
       7           0 :                 field(field), extends(extends), origin(0, 0, 0), reflective(false) {
       8             : 
       9           0 : }
      10             : 
      11           1 : PeriodicMagneticField::PeriodicMagneticField(ref_ptr<MagneticField> field,
      12           1 :                 const Vector3d &extends, const Vector3d &origin, bool reflective) :
      13           1 :                 field(field), extends(extends), origin(origin), reflective(reflective) {
      14             : 
      15           1 : }
      16             : 
      17           0 : Vector3d &PeriodicMagneticField::getOrigin() {
      18           0 :         return origin;
      19             : }
      20             : 
      21           0 : void PeriodicMagneticField::setOrigin(const Vector3d &origin) {
      22             :         this->origin = origin;
      23           0 : }
      24             : 
      25           0 : Vector3d &PeriodicMagneticField::getExtends() {
      26           0 :         return extends;
      27             : }
      28             : 
      29           0 : void PeriodicMagneticField::setExtends(const Vector3d &origin) {
      30             :         this->extends = extends;
      31           0 : }
      32             : 
      33           0 : bool PeriodicMagneticField::isReflective() {
      34           0 :         return reflective;
      35             : }
      36             : 
      37           0 : void PeriodicMagneticField::setReflective(bool reflective) {
      38           0 :         this->reflective = reflective;
      39           0 : }
      40             : 
      41           3 : Vector3d PeriodicMagneticField::getField(const Vector3d &position) const {
      42           3 :         Vector3d n = ((position - origin) / extends).floor();
      43             :         Vector3d p = position - origin - n * extends;
      44             : 
      45           3 :         if (reflective) {
      46           3 :                 long mx = (long) ::fabs(n.x) % 2;
      47           3 :                 if (mx == 1)
      48           2 :                         p.x = extends.x - p.x;
      49           3 :                 long my = (long) ::fabs(n.y) % 2;
      50           3 :                 if (my == 1)
      51           0 :                         p.y = extends.y - p.y;
      52           3 :                 long mz = (long) ::fabs(n.z) % 2;
      53           3 :                 if (mz == 1)
      54           2 :                         p.z = extends.z - p.z;
      55             :         }
      56             : 
      57           3 :         return field->getField(p);
      58             : }
      59             : 
      60           3 : void MagneticFieldList::addField(ref_ptr<MagneticField> field) {
      61           3 :         fields.push_back(field);
      62           3 : }
      63             : 
      64           1 : Vector3d MagneticFieldList::getField(const Vector3d &position) const {
      65             :         Vector3d b;
      66           4 :         for (int i = 0; i < fields.size(); i++)
      67           3 :                 b += fields[i]->getField(position);
      68           1 :         return b;
      69             : }
      70             : 
      71           1 : MagneticFieldEvolution::MagneticFieldEvolution(ref_ptr<MagneticField> field,
      72           1 :         double m) :
      73           1 :         field(field), m(m) {
      74           1 : }
      75             : 
      76           2 : Vector3d MagneticFieldEvolution::getField(const Vector3d &position,
      77             :         double z) const {
      78           2 :         return field->getField(position) * pow(1+z, m);
      79             : }
      80             : 
      81           2 : Vector3d MagneticDipoleField::getField(const Vector3d &position) const {
      82             :                 Vector3d r = (position - origin);
      83           2 :                 Vector3d unit_r = r.getUnitVector();
      84             :                 
      85           2 :                 if (r.getR() == 0) { // singularity
      86             :                         return moment * 2 * mu0 / 3;
      87             :                 }
      88           2 :                 return (unit_r * (unit_r.dot(moment)) * 3 - moment) / pow(r.getR() / radius, 3) * mu0 / (4*M_PI);
      89             : }
      90             : 
      91             : #ifdef CRPROPA_HAVE_MUPARSER
      92           1 : RenormalizeMagneticField::RenormalizeMagneticField(ref_ptr<MagneticField> field,
      93           1 :                 std::string expression) :
      94           2 :                 field(field), expression(expression) {
      95             : 
      96           1 :         p =  new mu::Parser();
      97           1 :         p->DefineVar("B", &Bmag);
      98           1 :         p->DefineConst("tesla", tesla);
      99           1 :         p->DefineConst("gauss", gauss);
     100           1 :         p->DefineConst("muG", muG);
     101           1 :         p->DefineConst("nG", nG);
     102           1 :         p->SetExpr(expression);
     103           1 : }
     104             : 
     105           1 : Vector3d RenormalizeMagneticField::getField(const Vector3d &position) {
     106           1 :         Vector3d B = field->getField(position);
     107           1 :         Bmag = B.getR();
     108           1 :         return B * p->Eval();
     109             : }
     110             : #endif
     111             : 
     112             : } // namespace crpropa