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 &extends) {
      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