Line data    Source code

       1             : #include "crpropa/massDistribution/Ferriere.h"
       2             : #include "crpropa/Common.h"
       3             : 
       4             : #include "kiss/logger.h"
       5             : 
       6             : #include <sstream>
       7             : 
       8             : namespace crpropa {
       9             : 
      10          13 : Vector3d Ferriere::CMZTransformation(const Vector3d &position) const {
      11             :         // set galactocentric coordinate system with the Sun at (-8.5,0.,0.) instead of (8.5, 0, 0) to be consistand with JF12 implementation
      12          13 :         double x = -position.x;
      13          13 :         double y = -position.y;
      14             : 
      15             :         double xC = -50*pc;             //offset
      16             :         double yC = 50*pc;
      17             :         double sinTc = sin(70.*deg);
      18             :         double cosTc = cos(70.*deg);
      19             : 
      20             :         Vector3d pos;
      21          13 :         pos.x = (x - xC)*cosTc + (y - yC)*sinTc;
      22          13 :         pos.y = -(x - xC)*sinTc + (y - yC)*cosTc;
      23          13 :         pos.z = position.z;
      24             : 
      25          13 :         return pos;
      26             : }
      27             : 
      28          13 : Vector3d Ferriere::DiskTransformation(const Vector3d &position) const {
      29             :         // set galactocentric coordinate system with the Sun at (-8.5,0.,0.) instead of (8.5, 0, 0) to be consistand with JF12 implementation
      30          13 :         double x = -position.x;
      31          13 :         double y = - position.y;
      32          13 :         double z = position.z;
      33             : 
      34             :         double alphaD = 13.5*deg;  // rotation arround x-axis
      35             :         double sinAd = sin(alphaD);
      36             :         double cosAd = cos(alphaD);
      37             :         double betaD = 20.*deg;  // rotation arround y'-axis
      38             :         double sinBd = sin(betaD);
      39             :         double cosBd = cos(betaD);
      40             :         double TettaD = 48.5*deg;  // rotation arround x"-axis
      41             :         double sinTd = sin(TettaD);
      42             :         double cosTd = cos(TettaD);
      43             : 
      44             :         Vector3d pos;
      45             : 
      46          13 :         pos.x = x*cosBd*cosTd - y*(sinAd*sinBd*cosTd -cosAd*sinTd)-z*(cosAd*sinBd*cosTd +sinAd*sinTd);
      47             : 
      48          13 :         pos.y =  -x*cosBd*sinTd;
      49          13 :         pos.y += y*(sinAd*sinBd*sinTd +cosAd*cosTd);
      50          13 :         pos.y += z*(cosAd*sinBd*sinTd -sinAd*cosTd);
      51             : 
      52          13 :         pos.z = x*sinBd;
      53          13 :         pos.z += y*sinAd*cosBd;
      54          13 :         pos.z += z*cosAd*cosBd;
      55             : 
      56          13 :         return pos;
      57             : }
      58             : 
      59          12 : double Ferriere::getHIDensity(const Vector3d &position) const {
      60             :         double n = 0;
      61          12 :         double R = sqrt(position.x*position.x+position.y*position.y);
      62             : 
      63          12 :         if(R<3*kpc)
      64             :         {
      65             :                 // density at center
      66           6 :                 Vector3d pos = CMZTransformation(position);  // coordinate trafo
      67           6 :                 double x = pos.x/pc;  // all units in pc
      68           6 :                 double y = pos.y/pc;
      69           6 :                 double z = pos.z/pc;
      70             : 
      71           6 :                 double A = sqrt(x*x+2.5*2.5*y*y);
      72           6 :                 double nCMZ = 8.8/ccm*exp(-pow_integer<4>((A-125.)/137))*exp(-pow_integer<2>(z/54.));
      73             : 
      74             :                 // density in disk
      75           6 :                 pos = DiskTransformation(position);  // Corrdinate Trafo
      76           6 :                 x = pos.x/pc;  // all units in pc
      77           6 :                 y = pos.y/pc;
      78           6 :                 z = pos.z/pc;
      79             : 
      80           6 :                 A = sqrt(x*x+3.1*3.1*y*y);
      81           6 :                 double nDisk = 0.34/ccm*exp(-pow_integer<4>((A-1200.)/438.))*exp(-pow_integer<2>(z/120));
      82             : 
      83           6 :                 n = nCMZ + nDisk;
      84             :         }
      85             :         else{  // outer region
      86           6 :                 double z = position.z/pc;
      87             :                 double a;
      88           6 :                 if(R<=Rsun){
      89             :                         a = 1;
      90             :                 } else {
      91           3 :                         a = R/Rsun;
      92             :                 }
      93             : 
      94           6 :                 double nCold = 0.859*exp(-pow_integer<2>(z/(127*a))); // cold HI
      95           6 :                 nCold += 0.047*exp(-pow_integer<2>(z/(318*a)));
      96           6 :                 nCold += 0.094*exp(-fabs(z)/(403*a));
      97           6 :                 nCold *= 0.340/ccm/(a*a);
      98             : 
      99           6 :                 double nWarm = (1.745 - 1.289/a)*exp(-pow_integer<2>(z/(127*a)));  // warm HI
     100           6 :                 nWarm += (0.473 - 0.070/a)*exp(-pow_integer<2>(z/(318*a)));
     101           6 :                 nWarm += (0.283 - 0.142/a)*exp(-fabs(z)/(403*a));
     102           6 :                 nWarm *= 0.226/ccm/a;
     103             : 
     104           6 :                 n = nWarm + nCold;
     105             :         }
     106             : 
     107          12 :         return n;
     108             : }
     109             : 
     110          12 : double Ferriere::getHIIDensity(const Vector3d &position) const {
     111             :         double n = 0;
     112          12 :         double R = sqrt(position.x*position.x+position.y*position.y);
     113             : 
     114          12 :         if(R< 3*kpc){   // inner
     115           6 :                 double x = position.x/pc;
     116           6 :                 double y = position.y/pc;
     117           6 :                 double z = position.z/pc;
     118             : 
     119             :                 // warm interstellar matter
     120           6 :                 double H = (x*x + pow_integer<2>(y+10.))/(145*145);
     121           6 :                 double nWIM = exp(-H)* exp(-pow_integer<2>(z+20.)/(26*26));
     122           6 :                 nWIM += 0.009*exp(-pow_integer<2>((R/pc-3700)/(0.5*3700)))*1/pow_integer<2>(cosh(z/140.));
     123           6 :                 nWIM += 0.005*cos(M_PI*R/pc*0.5/17000)*1/pow_integer<2>(cosh(z/950.));
     124           6 :                 nWIM *= 8.0/ccm;  // rescaling with density at center
     125             : 
     126             :                 //very hot interstellar matter
     127             :                 double alphaVH = 21.*deg;  // angle for very hot IM in radiant
     128             :                 double cosA = cos(alphaVH);
     129             :                 double sinA = sin(alphaVH);
     130           6 :                 double etta = y*cosA+z*sinA;  // coordinate transformation for VHIM along major axis
     131           6 :                 double chi = -y*sinA+z*cosA;
     132             : 
     133           6 :                 double nVHIM = 0.29/ccm*exp(-((x*x+etta*etta)/(162.*162.)+chi*chi/(90*90)));
     134             : 
     135           6 :                 n = nWIM + nVHIM;
     136             :         } else {  // outer region
     137           6 :                 double z = position.z;
     138             : 
     139           6 :                 double nWarm = 0.0237/ccm*exp(-(R*R-Rsun*Rsun)/pow_integer<2>(37*kpc))*exp(-fabs(z)/(1*kpc));
     140           6 :                 nWarm += 0.0013/ccm*exp(-(pow_integer<2>(R-4*kpc)-pow_integer<2>(Rsun-4*kpc))/pow_integer<2>(2*kpc))*exp(-fabs(z)/(150*pc));
     141             : 
     142           6 :                 double nHot = 0.12*exp(-(R-Rsun)/(4.9*kpc));
     143           6 :                 nHot += 0.88*exp(-(pow_integer<2>(R-4.5*kpc)-pow_integer<2>(Rsun-4.5*kpc))/pow_integer<2>(2.9*kpc));
     144           6 :                 nHot *= pow(R/Rsun, -1.65);
     145           6 :                 nHot *= exp(-fabs(z)/((1500*pc)*pow(R/Rsun,1.65)));
     146           6 :                 nHot *= 4.8e-4/ccm;
     147             : 
     148           6 :                 n= nWarm + nHot;
     149             :         }
     150          12 :         return n;
     151             : }
     152             : 
     153          12 : double Ferriere::getH2Density(const Vector3d &position) const{
     154             :         double n=0;
     155          12 :         double R=sqrt(position.x*position.x+position.y*position.y);
     156             : 
     157          12 :         if(R<3*kpc) {
     158             :                 // density at center
     159           6 :                 Vector3d pos =CMZTransformation(position);  // coord trafo for CMZ
     160           6 :                 double x = pos.x/pc;  // all units in pc
     161           6 :                 double y = pos.y/pc;
     162           6 :                 double z = pos.z/pc;
     163             : 
     164           6 :                 double A = sqrt(x*x+pow(2.5*y,2));  // ellipticity
     165           6 :                 double nCMZ = exp(-pow((A-125.)/137.,4))*exp(-pow(z/18.,2));
     166           6 :                 nCMZ *= 150/ccm;  // rescaling
     167             : 
     168             :                 // density in disk
     169           6 :                 pos = DiskTransformation(position);  // coord trafo for disk
     170           6 :                 x=pos.x/pc;
     171           6 :                 y=pos.y/pc;
     172           6 :                 z=pos.z/pc;
     173             : 
     174           6 :                 A = sqrt(x*x+pow_integer<2>(3.1*y));
     175           6 :                 double nDISK = exp(-pow_integer<4>((A-1200)/438))*exp(-pow_integer<2>(z/42));
     176           6 :                 nDISK *= 4.8/ccm;  // rescaling
     177             : 
     178           6 :                 n = nCMZ + nDISK;
     179             :         } else {  // outer region
     180           6 :                 double z = position.z/pc;
     181           6 :                 n = pow(R/Rsun, -0.58);
     182           6 :                 n *= exp(-(pow_integer<2>(R-4.5*kpc)-pow_integer<2>(Rsun-4.5*kpc))/pow_integer<2>(2.9*kpc));
     183           6 :                 n *= exp(-pow_integer<2>(z/(81*pow(R/Rsun,0.58))));
     184           6 :                 n *= 0.58/ccm;  // rescaling
     185             :         }
     186             : 
     187          12 :         return n;
     188             : }
     189             : 
     190           5 : double Ferriere::getDensity(const Vector3d &position) const{
     191             :         double n=0;
     192           5 :         if(isforHI){
     193           4 :                 n += getHIDensity(position);
     194             :         }
     195           5 :         if(isforHII){
     196           4 :                 n+=getHIIDensity(position);
     197             :         }
     198           5 :         if(isforH2){
     199           4 :                 n+=getH2Density(position);
     200             :         }
     201             :         // check if all densities are deactivated and raise warning if so
     202           5 :         if((isforHI || isforHII || isforH2) == false){
     203           2 :                 KISS_LOG_WARNING
     204           1 :                         << "\nCalled getDensity on fully deactivated Ferriere \n"
     205           1 :                         << "gas density model. The total density is set to 0.";
     206             :         }
     207             : 
     208           5 :         return n;
     209             : }
     210             : 
     211           5 : double Ferriere::getNucleonDensity(const Vector3d &position) const{
     212             :         double n=0;
     213           5 :         if(isforHI){
     214           4 :                 n += getHIDensity(position);
     215             :         }
     216           5 :         if(isforHII){
     217           4 :                 n+=getHIIDensity(position);
     218             :         }
     219           5 :         if(isforH2){
     220           4 :                 n+= 2*getH2Density(position);
     221             :         }
     222             : 
     223             :         // check if all densities are deactivated and raise warning if so
     224           5 :         if((isforHI || isforHII || isforH2) == false){
     225           2 :                 KISS_LOG_WARNING
     226           1 :                         << "\nCalled getNucleonDensity on fully deactivated Ferriere \n"
     227           1 :                         << "gas density model. The total density is set to 0.";
     228             :         }
     229             : 
     230           5 :         return n;
     231             : }
     232             : 
     233           1 : void Ferriere::setIsForHI(bool HI){
     234           1 :         isforHI = HI;
     235           1 : }
     236             : 
     237           1 : void Ferriere::setIsForHII(bool HII){
     238           1 :         isforHII = HII;
     239           1 : }
     240             : 
     241           1 : void Ferriere::setIsForH2(bool H2){
     242           1 :         isforH2 = H2;
     243           1 : }
     244             : 
     245           4 : bool Ferriere::getIsForHI(){
     246           4 :         return isforHI;
     247             : }
     248             : 
     249           4 : bool Ferriere::getIsForHII(){
     250           4 :         return isforHII;
     251             : }
     252             : 
     253           4 : bool Ferriere::getIsForH2(){
     254           4 :         return isforH2;
     255             : }
     256             : 
     257           0 : std::string Ferriere::getDescription() {
     258           0 :         std::stringstream s;
     259           0 :         s << "Density model Ferriere 2007:\n";
     260           0 :         s<< "HI component is ";
     261           0 :         if(!isforHI)
     262           0 :                 s<< "not ";
     263           0 :         s<< "active.\nHII component is ";
     264           0 :         if(!isforHII)
     265           0 :                 s<< "not ";
     266           0 :         s<<"active.\nH2 component is ";
     267           0 :         if(!isforH2)
     268           0 :                 s<<"not ";
     269           0 :         s<<"active.";
     270           0 :         return s.str();
     271           0 : }
     272             : 
     273             : }  // namespace crpropa