Line data    Source code

       1              : #include "crpropa/massDistribution/Nakanishi.h"
       2              : #include "crpropa/Common.h"
       3              : 
       4              : #include "kiss/logger.h"
       5              : 
       6              : #include <sstream>
       7              : 
       8              : namespace crpropa {
       9              : 
      10            8 : double Nakanishi::getHIScaleheight(const Vector3d &position) const {
      11            8 :         double R = sqrt(pow_integer<2>(position.x)+pow_integer<2>(position.y));      // radius in galactic plane
      12            8 :         double scaleheight = 1.06*pc*(116.3 +19.3*R/kpc+4.1*pow_integer<2>(R/kpc)-0.05*pow_integer<3>(R/kpc));
      13            8 :         return scaleheight;
      14              :         }
      15              : 
      16            8 : double Nakanishi::getHIPlanedensity(const Vector3d &position) const {
      17            8 :         double R = sqrt(pow_integer<2>(position.x)+pow_integer<2>(position.y));      // radius in galactic plane
      18            8 :         double planedensity = 0.94/ccm*(0.6*exp(-R/(2.4*kpc))+0.24*exp(-pow_integer<2>((R-9.5*kpc)/(4.8*kpc))));
      19            8 :         return planedensity;
      20              :         }
      21              : 
      22              : 
      23            8 : double Nakanishi::getH2Scaleheight(const Vector3d &position) const {
      24            8 :         double R = sqrt(pow_integer<2>(position.x)+ pow_integer<2>(position.y));  // radius in galactic plane
      25            8 :         double scaleheight = 1.06*pc*( 10.8*exp(0.28*R/kpc)+42.78);
      26            8 :         return scaleheight;
      27              : }
      28              : 
      29            8 : double Nakanishi::getH2Planedensity(const Vector3d &position) const {
      30            8 :         double R = sqrt(pow_integer<2>(position.x)+pow_integer<2>(position.y));  // radius in galactic plane
      31            8 :         double planedensity =0.94/ccm*(11.2*exp(-R*R/(0.874*kpc*kpc)) +0.83*exp(-pow_integer<2>((R-4*kpc)/(3.2*kpc))));
      32            8 :         return planedensity;
      33              : }
      34              : 
      35            6 : double Nakanishi::getHIDensity(const Vector3d &position) const {
      36              :         double n = 0;  // density
      37            6 :         double planedensity = getHIPlanedensity(position);
      38            6 :         double scaleheight = getHIScaleheight(position);
      39            6 :         n = planedensity*pow(0.5,pow_integer<2>(position.z/scaleheight));
      40              : 
      41            6 :         return n;
      42              : }
      43              : 
      44            6 : double Nakanishi::getH2Density(const Vector3d &position) const {
      45              :         double n = 0;  // density
      46            6 :         double planedensity = getH2Planedensity(position);
      47            6 :         double scaleheight = getH2Scaleheight(position);
      48            6 :         n = planedensity*pow(0.5,pow_integer<2>(position.z/scaleheight));
      49              : 
      50            6 :         return n;
      51              : }
      52              : 
      53            3 : double Nakanishi::getDensity(const Vector3d &position) const {
      54              :         double n = 0;
      55            3 :         if(isforHI)
      56            2 :                 n += getHIDensity(position);
      57            3 :         if(isforH2)
      58            2 :                 n += getH2Density(position);
      59              :         
      60              :         // check if all densities are deactivated and raise warning if so
      61            3 :         if((isforHI || isforH2) == false){
      62            2 :                 KISS_LOG_WARNING
      63            1 :                         << "\n"<<"Called getDensity on fully deactivated Nakanishi \n"
      64            1 :                         << "gas density model. The total density is set to 0.";
      65              :         }
      66            3 :         return n;
      67              : }
      68              : 
      69            3 : double Nakanishi::getNucleonDensity(const Vector3d &position) const {
      70              :         double n = 0;
      71            3 :         if(isforHI)
      72            2 :                 n += getHIDensity(position);
      73            3 :         if(isforH2)
      74            2 :                 n += 2*getH2Density(position);  // weight 2 for molecular hydrogen
      75              : 
      76              :         // check if all densities are deactivated and raise warning if so
      77            3 :         if((isforHI || isforH2) == false){
      78            2 :                 KISS_LOG_WARNING
      79            1 :                         << "\n"<<"Called getNucleonDensity on fully deactivated Nakanishi \n"
      80            1 :                         << "gas density model. The total density is set to 0.";
      81              :         }
      82              : 
      83            3 :         return n;
      84              : }
      85              : 
      86            3 : bool Nakanishi::getIsForHI() {
      87            3 :         return isforHI;
      88              : }
      89              : 
      90            1 : bool Nakanishi::getIsForHII() {
      91            1 :         return isforHII;
      92              : }
      93            3 : bool Nakanishi::getIsForH2() {
      94            3 :         return isforH2;
      95              : }
      96              : 
      97            1 : void Nakanishi::setIsForHI(bool HI) {
      98            1 :         isforHI = HI;
      99            1 : }
     100              : 
     101            1 : void Nakanishi::setIsForH2(bool H2) {
     102            1 :         isforH2 = H2;
     103            1 : }
     104              : 
     105            0 : std::string Nakanishi::getDescription() {
     106            0 :         std::stringstream s;
     107            0 :         s << "Density model Nakanishi:\n";
     108            0 :         s<< "HI component is ";
     109            0 :         if(isforHI==false)
     110            0 :                 s<< "not ";
     111            0 :         s<< "active.\nH2 component is ";
     112            0 :         if(isforH2==false)
     113            0 :                 s<<"not ";
     114            0 :         s<<"active.\nNakanishi has no HII component.";
     115              : 
     116            0 :         return s.str();
     117            0 : }
     118              : 
     119              : }  // namespace crpropa