Line data    Source code

       1              : #include "crpropa/massDistribution/Massdistribution.h"
       2              : #include <sstream>
       3              : namespace crpropa {
       4              : 
       5            2 : void DensityList::addDensity(ref_ptr<Density> dens) {
       6            2 :         DensityList.push_back(dens);
       7            2 : }
       8              : 
       9            1 : double DensityList::getDensity(const Vector3d &position) const {
      10              :         double n = 0.;
      11            3 :         for (int i = 0; i < DensityList.size(); i++)
      12            2 :                 n += DensityList[i]->getDensity(position);
      13            1 :         return n;
      14              : }
      15              : 
      16            1 : double DensityList::getHIDensity(const Vector3d &position) const {
      17              :         double n = 0.;
      18            3 :         for (int i = 0; i < DensityList.size(); i++)
      19            2 :                 n += DensityList[i]->getHIDensity(position);
      20            1 :         return n;
      21              : }
      22              : 
      23            1 : double DensityList::getHIIDensity(const Vector3d &position) const {
      24              :         double n = 0.;
      25            3 :         for (int i = 0; i < DensityList.size(); i++)
      26            2 :                 n += DensityList[i]->getHIIDensity(position);
      27            1 :         return n;
      28              : }
      29              : 
      30            1 : double DensityList::getH2Density(const Vector3d &position) const {
      31              :         double n = 0.;
      32            3 :         for (int i = 0; i < DensityList.size(); i++)
      33            2 :                 n += DensityList[i]->getH2Density(position);
      34            1 :         return n;
      35              : }
      36              : 
      37            1 : double DensityList::getNucleonDensity(const Vector3d &position) const {
      38              :         double n = 0.;
      39            3 :         for (int i = 0; i < DensityList.size(); i++)
      40            2 :                 n += DensityList[i]->getNucleonDensity(position);
      41            1 :         return n;
      42              : }
      43              : 
      44            0 : std::string DensityList::getDescription() {
      45            0 :         std::stringstream ss; 
      46            0 :         ss << "DensityList with " << DensityList.size() << " modules: \n";
      47            0 :         for (int i = 0; i < DensityList.size(); i++) {
      48            0 :                 ss << "density " << i + 1 << ": " << DensityList[i] -> getDescription();
      49              :         }
      50              :         
      51            0 :         return ss.str();
      52            0 : }
      53              : 
      54              : // ----------- DensityGrid -----------------------------------------------------------------
      55              : 
      56            2 : DensityGrid::DensityGrid(ref_ptr<Grid1f> grid, bool isForHI, bool isForHII, bool isForH2) : 
      57            2 :         grid(grid), isForHI(isForHI), isForHII(isForHII), isForH2(isForH2) {
      58            2 :                 checkAndWarn();
      59            2 :         }
      60              : 
      61            5 : void DensityGrid::checkAndWarn() {
      62            5 :         bool allDeactivated = (isForHI == false) && (isForHII == false) && (isForH2 == false);
      63              :         if (allDeactivated) {
      64            0 :                 KISS_LOG_WARNING << "DensityGrid has all types deactivated."
      65            0 :                         << "In this case all output will be n = 0. \n"
      66            0 :                         << "Please activate the intended particle type. \n";
      67              :         }
      68            5 : }
      69              : 
      70            3 : double DensityGrid::getHIDensity(const Vector3d &position) const {
      71            3 :         if (isForHI)
      72            3 :                 return grid -> interpolate(position);
      73              :         else 
      74              :                 return 0.;
      75              : }
      76              : 
      77            3 : double DensityGrid::getHIIDensity(const Vector3d &position) const {
      78            3 :         if (isForHII) 
      79            0 :                 return grid -> interpolate(position);
      80              :         else
      81              :                 return 0.;
      82              : }
      83              : 
      84            3 : double DensityGrid::getH2Density(const Vector3d &position) const {
      85            3 :         if (isForH2)
      86            0 :                 return grid -> interpolate(position);
      87              :         else
      88              :                 return 0.;
      89              : }
      90              : 
      91            1 : double DensityGrid::getDensity(const Vector3d &position) const {
      92              :         double n = 0;
      93            1 :         n += getHIDensity(position);
      94            1 :         n += getHIIDensity(position);
      95            1 :         n += getH2Density(position);
      96              : 
      97            1 :         return n;
      98              : }
      99              : 
     100            1 : double DensityGrid::getNucleonDensity(const Vector3d &position) const {
     101              :         double n = 0;
     102            1 :         n += getHIDensity(position);
     103            1 :         n += getHIIDensity(position);
     104            1 :         n += 2 * getH2Density(position);
     105              : 
     106            1 :         return n;
     107              : }
     108              : 
     109            2 : bool DensityGrid::getIsForHI() {
     110            2 :         return isForHI;
     111              : }
     112              : 
     113            2 : bool DensityGrid::getIsForHII() {
     114            2 :         return isForHII;
     115              : }
     116              : 
     117            2 : bool DensityGrid::getIsForH2() {
     118            2 :         return isForH2;
     119              : }
     120              : 
     121            1 : void DensityGrid::setIsForHI(bool b) {
     122            1 :         isForHI = b;
     123            1 :         checkAndWarn();
     124            1 : }
     125              : 
     126            1 : void DensityGrid::setIsForHII(bool b) {
     127            1 :         isForHII = b;
     128            1 :         checkAndWarn();
     129            1 : }
     130              : 
     131            1 : void DensityGrid::setIsForH2(bool b) {
     132            1 :         isForH2 = b;
     133            1 :         checkAndWarn();
     134            1 : }
     135              : 
     136            0 : void DensityGrid::setGrid(ref_ptr<Grid1f> grid) {
     137            0 :         this->grid = grid;
     138            0 : }
     139              : 
     140            0 : std::string DensityGrid::getDescription() {
     141            0 :         std::stringstream ss;
     142            0 :         ss << "Density in a given grid \n"; 
     143            0 :         return ss.str();
     144            0 : }
     145              : 
     146              : } //namespace crpropa