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