Commit dea2661a722bcd820e093004c32650e59c130dee

Authored by Jürgen Knödlseder
1 parent e156834b

Add HealPix flux computation test.

Showing 1 changed file with 62 additions and 0 deletions
test/test_healpix_flux.py 0 → 100755
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script tests the accuracy of the HealPix flux computation. It sets
  4 +# a single pixel of the HealPix map, projects the HealPix map into a CAR
  5 +# map to compute the flux in the pixel in an indpendent way, and compare
  6 +# this flux estimates to that one derived from the flux() method.
  7 +#
  8 +# Note that this test relies on the correctness of the CAR flux()
  9 +# computation.
  10 +#
  11 +# --------------------------------------------------------------------------
  12 +#
  13 +# Copyright (C) 2015 Juergen Knoedlseder
  14 +#
  15 +# This program is free software: you can redistribute it and/or modify
  16 +# it under the terms of the GNU General Public License as published by
  17 +# the Free Software Foundation, either version 3 of the License, or
  18 +# (at your option) any later version.
  19 +#
  20 +# This program is distributed in the hope that it will be useful,
  21 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  23 +# GNU General Public License for more details.
  24 +#
  25 +# You should have received a copy of the GNU General Public License
  26 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  27 +#
  28 +# ==========================================================================
  29 +import gammalib
  30 +
  31 +# Allocate HealPix map
  32 +nside = 16
  33 +npix = 12 * nside*nside
  34 +print("Pixels %d" % npix)
  35 +
  36 +# Loop over all pixels
  37 +for index in range(npix):
  38 +
  39 + # Allocate HealPix map
  40 + map = gammalib.GSkymap("GAL", nside, "RING")
  41 +
  42 + # Create WCS map, set Healpix index and merge into WCS map
  43 + wcs = gammalib.GSkymap("CAR", "GAL", 0.0, 0.0, -0.1, 0.1, 3600, 1800)
  44 + map[index] = 1.0
  45 + wcs += map
  46 +
  47 + # Get total flux in WCS map (in HealPix pixel only)
  48 + total = 0.0
  49 + for i in range(wcs.npix()):
  50 + if map.dir2inx(wcs.inx2dir(i)) == index:
  51 + total += wcs.flux(i)
  52 + else:
  53 + wcs[i] = 0
  54 +
  55 + # Save map
  56 + wcs.save("healpix.fits", True)
  57 +
  58 + # Print results
  59 + print("Index %d" % index)
  60 + print("Flux in pixel from integration .....: %f" % total)
  61 + print("Flux in pixel from flux() method ...: %f (%f%%)" % \
  62 + (map.flux(index),(map.flux(index)-total)/total*100.0))
0 63 \ No newline at end of file
... ...