Line data    Source code

       1             : //----------------------------------------------------------------------
       2             : // This file is part of PARSEC (http://physik.rwth-aachen.de/parsec)
       3             : // a parametrized simulation engine for cosmic rays.
       4             : //
       5             : // Copyright (C) 2011  Martin Erdmann, Peter Schiffer, Tobias Winchen
       6             : //                     RWTH Aachen University, Germany
       7             : // Contact: winchen@physik.rwth-aachen.de
       8             : //
       9             : //  This program is free software: you can redistribute it and/or
      10             : //  modify it under the terms of the GNU General Public License as
      11             : //  published by the Free Software Foundation, either version 3 of
      12             : //  the License, or (at your option) any later version.
      13             : //
      14             : //  This program is distributed in the hope that it will be useful,
      15             : //  but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             : //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             : //  GNU General Public License for more details.
      18             : //
      19             : //  You should have received a copy of the GNU General Public License
      20             : //  along with this program. If not, see <http://www.gnu.org/licenses/>.
      21             : //----------------------------------------------------------------------
      22             : 
      23             : #include "crpropa/magneticLens/Pixelization.h"
      24             : #include "crpropa/Random.h"
      25             : 
      26             : namespace crpropa 
      27             : {
      28             : 
      29             :         healpix::T_Healpix_Base<healpix::int64> Pixelization::_healpix_nest = healpix::T_Healpix_Base<healpix::int64>(29, healpix::NEST);
      30             : 
      31             : 
      32           0 : uint8_t Pixelization::pix2Order(uint32_t pix)
      33             : {
      34           0 :         for (uint8_t i = 0; i < _nOrder_max; i++)
      35             :         {
      36           0 :                 if (pix == _nPix[i])
      37           0 :                         return i + 1;
      38             :         }
      39             :         return 0;
      40             : }
      41             : 
      42           0 : uint32_t Pixelization::nPix(uint8_t order)
      43             : {
      44           0 :         if (order > _nOrder_max)
      45             :         {
      46             :                 return 0;
      47             :         }
      48             :         else
      49             :         {
      50           0 :                 return _nPix[order - 1];
      51             :         }
      52             : }
      53             : 
      54       24585 : uint32_t Pixelization::direction2Pix(double longitude, double latitude) const
      55             : {
      56             :         healpix::vec3 v;
      57       24585 :         spherCo2Vec(longitude, latitude, v);
      58             :         try
      59             :         {
      60       24585 :                 uint32_t i = (uint32_t) _healpix->vec2pix(v);
      61       24585 :                 return i;
      62             :         }
      63           0 :         catch (healpix::PlanckError &e)
      64             :         {
      65           0 :                 std::cerr << "Healpix error triggered from direction2Pix(" << longitude << ", " << latitude  << ")\n";
      66           0 :                 std::cerr << " v = " << v.x <<", " << v.y << ", " <<  v.z << std::endl;
      67           0 :                 std::cerr << "\n The original exception reads:\n";
      68           0 :                 std::cerr << e.what() << std::endl;
      69           0 :                 throw;
      70           0 :         }
      71             : }
      72             : 
      73       36868 : void Pixelization::pix2Direction(uint32_t i, double &longitude,
      74             :                 double &latitude) const
      75             : {
      76             :         healpix::vec3 v;
      77             :         try{
      78       36868 :                 v = _healpix->pix2vec(i);
      79             :         }
      80           0 :         catch (healpix::PlanckError &e)
      81             :         {
      82           0 :                 std::cerr << "Healpix error triggered from pix2Direction(" << i << ", &longitude, &latitude " << ")\n";
      83           0 :                 std::cerr << "The original exception reads:\n";
      84           0 :                 std::cerr << e.what() << std::endl;
      85           0 :                 throw;
      86           0 :         }
      87             : 
      88       36868 :         vec2SphereCo(longitude, latitude, v);
      89       36868 : }
      90             : 
      91       24585 : void Pixelization::spherCo2Vec(double phi, double theta,
      92             :                 healpix::vec3 &V) const
      93             : {
      94       24585 :         V.x = cos(phi) * cos(theta);
      95       24585 :         V.y = sin(phi) * cos(theta);
      96       24585 :         V.z = sin(theta);
      97       24585 : }
      98             : 
      99       37289 : void Pixelization::vec2SphereCo(double &phi, double &theta,
     100             :                 const healpix::vec3 &V) const
     101             : {
     102       37289 :         theta = asin(V.z);
     103       37289 :         phi = healpix::safe_atan2(V.y, V.x);
     104       37289 : }
     105             : 
     106             : 
     107       49152 : double Pixelization::angularDistance(uint32_t i, uint32_t j) const
     108             : {
     109             :         healpix::vec3 v1, v2;
     110       49152 :         v1 = _healpix->pix2vec(i);
     111       49152 :         v2 = _healpix->pix2vec(j);
     112       49152 :         double s = v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
     113             :         // Failsafe for numerical inaccuracies
     114       49152 :         return ((s > 1) ? 0 : ((s < -1) ? M_PI : acos(s)));
     115             : }
     116             : 
     117         421 : void Pixelization::getRandomDirectionInPixel(uint32_t i, double &longitude, double &latitude) 
     118             : {
     119             :         
     120         421 :         uint64_t inest = _healpix->ring2nest(i);
     121         421 :         uint64_t nUp = 29 - _healpix->Order();
     122         421 :         uint64_t iUp = inest * pow(4, nUp);
     123         421 :         iUp += Random::instance().randInt64(pow(4, nUp));
     124             : 
     125         421 :         healpix::vec3 v = _healpix_nest.pix2vec(iUp);
     126             :         
     127         421 :         vec2SphereCo(longitude, latitude, v);
     128         421 : }
     129             : } // namespace