Line data    Source code

       1             : #ifndef CRPROPA_ADVECTIONFIELD_H
       2             : #define CRPROPA_ADVECTIONFIELD_H
       3             : 
       4             : 
       5             : #include <string>
       6             : #include <iostream>
       7             : #include <cmath>
       8             : #include <cstdlib>
       9             : #include <sstream>
      10             : 
      11             : #include "crpropa/Vector3.h"
      12             : #include "crpropa/Referenced.h"
      13             : #include "crpropa/Units.h"
      14             : 
      15             : namespace crpropa {
      16             : 
      17             : /**
      18             :  @class AdvectionField
      19             :  @brief Abstract base class for advection fields. These are used to model 
      20             :         the deterministic part of the Fokker-Planck equation. The getDivergence()
      21             :         method is used to model the adibatic cooling/heating.
      22             :  */
      23             : class AdvectionField: public Referenced {
      24             : public:
      25             :         virtual ~AdvectionField() {
      26             :         }
      27             :         virtual Vector3d getField(const Vector3d &position) const = 0;
      28             :         virtual double getDivergence(const Vector3d &position) const = 0;
      29             : };
      30             : 
      31             : 
      32             : /**
      33             :  @class AdvectionFieldList
      34             :  @brief Advection field decorator implementing a superposition of fields.
      35             :  */
      36           2 : class AdvectionFieldList: public AdvectionField {
      37             :         std::vector<ref_ptr<AdvectionField> > fields;
      38             : public:
      39             :         void addField(ref_ptr<AdvectionField> field);
      40             :         Vector3d getField(const Vector3d &position) const;
      41             :         double getDivergence(const Vector3d &position) const;
      42             : };
      43             : 
      44             : 
      45             : /**
      46             :  @class UniformAdvectionField
      47             :  @brief Advection field with one velocity/advection-field vector.
      48             :  */
      49           1 : class UniformAdvectionField: public AdvectionField {
      50             :         Vector3d value;
      51             : public:
      52             :         UniformAdvectionField(const Vector3d &value);
      53             :         Vector3d getField(const Vector3d &position) const;
      54             :         double getDivergence(const Vector3d &position) const;
      55             : 
      56             :         std::string getDescription() const;
      57             : };
      58             : 
      59             : 
      60             : /**
      61             : @class ConstantSphericalAdvectionField
      62             : @brief Spherical advection field with a constant wind speed
      63             : */
      64             : 
      65           1 : class ConstantSphericalAdvectionField: public AdvectionField {
      66             :         Vector3d origin; //origin of the advection sphere
      67             :         double vWind; // wind velocity
      68             : public:
      69             :         /** Constructor
      70             :          @param origin  Origin of the advection field
      71             :          @param vWind   Constant wind velocity
      72             : 
      73             : */
      74             : 
      75             :         ConstantSphericalAdvectionField(const Vector3d origin, double vWind);
      76             :         Vector3d getField(const Vector3d &position) const;
      77             :         double getDivergence(const Vector3d &position) const;
      78             : 
      79             :         void setOrigin(const Vector3d origin);
      80             :         void setVWind(double vMax);
      81             : 
      82             :         Vector3d getOrigin() const;
      83             :         double getVWind() const;
      84             : 
      85             :         std::string getDescription() const;
      86             : 
      87             :         
      88             : };
      89             :         
      90             : /**
      91             :  @class SphericalAdvectionField
      92             :  @brief Spherical advection with a exponentially increasing and
      93             :         exponentially constant velocity.
      94             : */
      95             : 
      96           1 : class SphericalAdvectionField: public AdvectionField {
      97             :         Vector3d origin; //origin of the advection sphere
      98             :         double radius; //radius of the advection sphere
      99             :         double vMax; // maximum wind velocity
     100             :         double tau; // transition distance
     101             :         double alpha; //tuning parameter
     102             : public:
     103             :         /** Constructor
     104             :         @param origin   Origin of the advection sphere
     105             :         @param radius   Radius of the advection sphere
     106             :         @param vMax     Maximum wind velocity
     107             :         @param tau      Transition distance
     108             :         @param alpha    Tuning parameter
     109             : */
     110             :         SphericalAdvectionField(const Vector3d origin, double radius, double vMax, double tau, double alpha);
     111             :         Vector3d getField(const Vector3d &position) const;
     112             :         double getDivergence(const Vector3d &position) const;
     113             : 
     114             :         double getV(const double &r) const;
     115             : 
     116             :         void setOrigin(const Vector3d origin);
     117             :         void setRadius(double radius);
     118             :         void setVMax(double vMax);
     119             :         void setTau(double tau);
     120             :         void setAlpha(double alpha);
     121             : 
     122             :         Vector3d getOrigin() const;
     123             :         double getRadius() const;
     124             :         double getVMax() const;
     125             :         double getTau() const;
     126             :         double getAlpha() const;
     127             :         
     128             :         std::string getDescription() const;
     129             : };
     130             : 
     131             : /**
     132             :  @class OneDimensionalCartesianShock
     133             :  @brief Advection field in x-direction with shock at x = 0 and width lShock approximated by tanh() 
     134             :                 with variable compression ratio vUp/vDown
     135             :  */
     136             : class OneDimensionalCartesianShock: public AdvectionField {
     137             :         double compressionRatio; //compression ratio of shock
     138             :         double vUp; //upstream velocity 
     139             :         double lShock; //shock width
     140             : public:
     141             : /** Constructor
     142             :         @param compressionRatio //compression ratio of shock
     143             :         @param vUp //upstream velocity 
     144             :         @param lShock //shock width
     145             : */
     146             :         OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock);
     147             :         Vector3d getField(const Vector3d &position) const;
     148             :         double getDivergence(const Vector3d &position) const;
     149             : 
     150             :         void setComp(double compressionRatio);
     151             :         void setVup(double vUp);
     152             :         void setShockwidth(double lShock);
     153             : 
     154             :         double getComp() const;
     155             :         double getVup() const; 
     156             :         double getShockwidth() const;
     157             : 
     158             :         std::string getDescription() const;
     159             : };
     160             : 
     161             : /**
     162             :  @class OneDimensionalSphericalShock
     163             :  @brief Advection field in x-direction with shock at rShock and width lShock approximated by tanh() 
     164             :                 with variable compression ratio ratio vUp/vDown
     165             :  */
     166             : class OneDimensionalSphericalShock: public AdvectionField {
     167             :         double compressionRatio;        //compression ratio of shock
     168             :         double vUp;     //upstream velocity 
     169             :         double lShock;  //shock width
     170             :         double rShock;  //shock radius
     171             :         bool coolUpstream;      //flag for upstream cooling
     172             : public:
     173             : /** Constructor
     174             :         @param compressionRatio //compression ratio of shock
     175             :         @param vUp      //upstream velocity 
     176             :         @param lShock   //shock width
     177             :         @param rShock   //shock radius
     178             :         @param coolUpstream //flag for upstream cooling
     179             : */
     180             :         OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream);
     181             :         Vector3d getField(const Vector3d &position) const;
     182             :         double getDivergence(const Vector3d &position) const;
     183             : 
     184             :         void setComp(double compressionRatio);
     185             :         void setVup(double vUp);
     186             :         void setShockwidth(double lShock);
     187             :         void setShockRadius(double rShock);
     188             :         void setCooling(bool coolUpstream);
     189             : 
     190             :         double getComp() const; 
     191             :         double getVup() const;
     192             :         double getShockwidth() const;
     193             :         double getShockRadius() const;
     194             :         bool getCooling() const;
     195             : 
     196             :         std::string getDescription() const;
     197             : };
     198             : 
     199             : /**
     200             :  @class ObliqueAdvectionShock
     201             :  @brief Advection field in x-y-direction with shock at x = 0 and width x_sh approximated by tanh() 
     202             :                 with variable compression ratio r_comp = vx_up/vx_down. The y component vy is not shocked 
     203             :                 and remains constant. 
     204             :  */
     205             : class ObliqueAdvectionShock: public AdvectionField {
     206             :         double compressionRatio; //compression ratio of shock
     207             :         double vXUp; //upstream velocity x-component
     208             :         double vY; //constant velocity y-component
     209             :         double lShock; //shock width
     210             :         
     211             : public:
     212             : /** Constructor
     213             :         @param compressionRatio //compression ratio of shock
     214             :         @param vXUp //upstream velocity x-component
     215             :         @param vY //constant velocity y-component
     216             :         @param lShock //shock width
     217             :         
     218             : */
     219             :         ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock);
     220             :         Vector3d getField(const Vector3d &position) const;
     221             :         double getDivergence(const Vector3d &position) const;
     222             : 
     223             :         void setComp(double compressionRatio);
     224             :         void setVup(double vXUp);
     225             :         void setVy(double vY);
     226             :         void setShockwidth(double lShock);
     227             : 
     228             :         double getComp() const; 
     229             :         double getVup() const;
     230             :         double getVy() const;
     231             :         double getShockwidth() const;
     232             :         
     233             :         std::string getDescription() const;
     234             : };
     235             : 
     236             : /**
     237             :  @class SphericalAdvectionShock
     238             :  @brief Spherical advection with a constant velocity for r<r_0
     239             :         at the the shock the velocity drops to v_0/4. followed by
     240             :         a decrease proportional to 1/r^2.
     241             : */
     242             : 
     243           1 : class SphericalAdvectionShock: public AdvectionField {
     244             :         Vector3d origin; // origin of the advection sphere
     245             :         double r_0; // position of the shock
     246             :         double v_0; // constant velocity
     247             :         double lambda; //transition width
     248             :         double r_rot; // normalization radius for rotation speed
     249             :         double v_phi; // rotation speed at r_rot
     250             : 
     251             : public:
     252             :         /** Constructor
     253             :         @param origin   Origin of the advection sphere
     254             :         @param r_0      Position of the shock
     255             :         @param v_0      Constant velocity (r<<r_o)
     256             :         @param lambda   Transition width / width of the shock
     257             : */
     258             :         SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double lambda);
     259             : 
     260             :         Vector3d getField(const Vector3d &position) const;
     261             :         double getDivergence(const Vector3d &position) const;
     262             : 
     263             :         double g(double R) const;
     264             :         double g_prime(double R) const;
     265             : 
     266             :         void setOrigin(const Vector3d Origin);
     267             :         void setR0(double r);
     268             :         void setV0(double v);
     269             :         void setLambda(double l);
     270             :         void setRRot(double r);
     271             :         void setAzimuthalSpeed(double vPhi);
     272             : 
     273             :         Vector3d getOrigin() const;
     274             :         double getR0() const;
     275             :         double getV0() const;
     276             :         double getLambda() const;
     277             :         /**
     278             :          * @param r Normalization radius for rotation speed
     279             :         */      
     280             :         double getRRot() const;
     281             :         /**
     282             :          * @param vPhi  Rotation speed at r_rot
     283             :         */      
     284             :         double getAzimuthalSpeed() const;
     285             : 
     286             :         std::string getDescription() const;
     287             : };
     288             : 
     289             : } // namespace crpropa
     290             : 
     291             : #endif // CRPROPA_ADVECTIONFIELD_H