Justine Devin / gammalib

Commit 69a208a2a73af0909bd4c162cca48f3bf252559b

Authored by Jürgen Knödlseder
1 parent 91bf6a00

Implement 8-wedge integration formula for WCS pixels

Inline Side-by-side
Showing 1 changed file with 90 additions and 56 deletions
src/sky/GSkymap.cpp
  View file @69a208a
... ... @@ -29,6 +29,7 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include "GException.hpp"
  32 +#include "GMath.hpp"
32 33 #include "GTools.hpp"
33 34 #include "GSkymap.hpp"
34 35 #include "GHealpix.hpp"
... ... @@ -851,13 +852,26 @@ double GSkymap::operator()(const GSkyDir&amp; dir, const int&amp; map) const
851 852 double x_size = std::abs(m_num_x * static_cast<GWcs*>(m_proj)->cdelt(0));
852 853 bool x_wrap = (x_size > 359.99);
853 854  
  855 + // Get pixel index that is to the left-top of the actual
  856 + // pixel. We take care of the special case of negative pixel
  857 + // indices which arise if we are in the first column or
  858 + // first row of the map.
  859 + int inx_x = int(pixel.x());
  860 + int inx_y = int(pixel.y());
  861 + if (pixel.x() < 0.0) {
  862 + inx_x--;
  863 + }
  864 + if (pixel.y() < 0.0) {
  865 + inx_y--;
  866 + }
  867 +
854 868 // Set left and right indices for interpolation. The left
855   - // index needs to non-negative and the right index needs
  869 + // index needs to be non-negative and the right index needs
856 870 // to be not larger that the number of pixels. We treat
857 871 // here also the special case of wrap around in the x
858 872 // axis that may occur if we have an allsky map.
859   - int inx_x_left = int(pixel.x());
860   - int inx_y_left = int(pixel.y());
  873 + int inx_x_left = inx_x;
  874 + int inx_y_left = inx_y;
861 875 int inx_x_right = inx_x_left + 1;
862 876 int inx_y_right = inx_y_left + 1;
863 877 if (inx_x_left < 0) {
... ... @@ -865,7 +879,7 @@ double GSkymap::operator()(const GSkyDir&amp; dir, const int&amp; map) const
865 879 inx_x_left += m_num_x;
866 880 }
867 881 else {
868   - inx_x_left = 0;
  882 + inx_x_left = 0;
869 883 }
870 884 }
871 885 if (inx_x_right >= m_num_x) {
... ... @@ -877,16 +891,16 @@ double GSkymap::operator()(const GSkyDir&amp; dir, const int&amp; map) const
877 891 }
878 892 }
879 893 if (inx_y_left < 0) {
880   - inx_y_left = 0;
  894 + inx_y_left = 0;
881 895 }
882 896 if (inx_y_right >= m_num_y) {
883 897 inx_y_right = m_num_y - 1;
884 898 }
885 899  
886 900 // Set weighting factors for interpolation
887   - double wgt_x_right = (pixel.x() - int(pixel.x()));
  901 + double wgt_x_right = (pixel.x() - inx_x);
888 902 double wgt_x_left = 1.0 - wgt_x_right;
889   - double wgt_y_right = (pixel.y() - int(pixel.y()));
  903 + double wgt_y_right = (pixel.y() - inx_y);
890 904 double wgt_y_left = 1.0 - wgt_y_right;
891 905  
892 906 // Compute skymap pixel indices for bi-linear interpolation
... ... @@ -1266,15 +1280,14 @@ double GSkymap::flux(const int&amp; index, const int&amp; map) const
1266 1280 *
1267 1281 * Returns the flux in the specified sky map @p pixel. The flux is computed
1268 1282 * by integrating the intensity over the solid angle subtended by the pixel.
  1283 + * Integration is done by dividing the pixel into wedges, computing the
  1284 + * intensitites at the cornes of each wedge, averaging these intensities,
  1285 + * and multiplying it with the solid angle of each wedge. This provides an
  1286 + * approximation of the true pixel flux which is accurate to better than
  1287 + * about 5%.
1269 1288 *
1270   - * For a HealPix pixelisation, the pixel is divided into 12 wedges and the
1271   - * flux is integrated in each wedge by computing the average of the corner
1272   - * intensities multiplied by the solid angle of the wedge. This leads to a
1273   - * precision of the order of 1-3% in the flux computation.
1274   - *
1275   - * For a WCS pixelisation, the pixel is divided into 4 quadrants and the flux
1276   - * is computed in each quandrant by multiplying the intensity in the centre
1277   - * of the quadrant with the solid angle of the quadrant.
  1289 + * For a HealPix pixelisation, the pixel is divided into 12 wedges. For a
  1290 + * WCS pixelisation, the pixel is divided into 8 wedges.
1278 1291 ***************************************************************************/
1279 1292 double GSkymap::flux(const GSkyPixel& pixel, const int& map) const
1280 1293 {
... ... @@ -1319,29 +1332,29 @@ double GSkymap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
1319 1332 double i12 = this->operator()(boundaries[11]);
1320 1333  
1321 1334 // Compute fluxes in wedges
1322   - double flux1 = (i0 + i1 + i2) / 3.0 *
  1335 + double flux1 = gammalib::onethird * (i0 + i1 + i2) *
1323 1336 solidangle(centre, boundaries[0], boundaries[4]);
1324   - double flux2 = (i0 + i2 + i3) / 3.0 *
  1337 + double flux2 = gammalib::onethird * (i0 + i2 + i3) *
1325 1338 solidangle(centre, boundaries[4], boundaries[8]);
1326   - double flux3 = (i0 + i3 + i4) / 3.0 *
  1339 + double flux3 = gammalib::onethird * (i0 + i3 + i4) *
1327 1340 solidangle(centre, boundaries[8], boundaries[1]);
1328   - double flux4 = (i0 + i4 + i5) / 3.0 *
  1341 + double flux4 = gammalib::onethird * (i0 + i4 + i5) *
1329 1342 solidangle(centre, boundaries[1], boundaries[5]);
1330   - double flux5 = (i0 + i5 + i6) / 3.0 *
  1343 + double flux5 = gammalib::onethird * (i0 + i5 + i6) *
1331 1344 solidangle(centre, boundaries[5], boundaries[9]);
1332   - double flux6 = (i0 + i6 + i7) / 3.0 *
  1345 + double flux6 = gammalib::onethird * (i0 + i6 + i7) *
1333 1346 solidangle(centre, boundaries[9], boundaries[2]);
1334   - double flux7 = (i0 + i7 + i8) / 3.0 *
  1347 + double flux7 = gammalib::onethird * (i0 + i7 + i8) *
1335 1348 solidangle(centre, boundaries[2], boundaries[6]);
1336   - double flux8 = (i0 + i8 + i9) / 3.0 *
  1349 + double flux8 = gammalib::onethird * (i0 + i8 + i9) *
1337 1350 solidangle(centre, boundaries[6], boundaries[10]);
1338   - double flux9 = (i0 + i9 + i10) / 3.0 *
  1351 + double flux9 = gammalib::onethird * (i0 + i9 + i10) *
1339 1352 solidangle(centre, boundaries[10], boundaries[3]);
1340   - double flux10 = (i0 + i10 + i11) / 3.0 *
  1353 + double flux10 = gammalib::onethird * (i0 + i10 + i11) *
1341 1354 solidangle(centre, boundaries[3], boundaries[7]);
1342   - double flux11 = (i0 + i11 + i12) / 3.0 *
  1355 + double flux11 = gammalib::onethird * (i0 + i11 + i12) *
1343 1356 solidangle(centre, boundaries[7], boundaries[11]);
1344   - double flux12 = (i0 + i12 + i1) / 3.0 *
  1357 + double flux12 = gammalib::onethird * (i0 + i12 + i1) *
1345 1358 solidangle(centre, boundaries[11], boundaries[0]);
1346 1359  
1347 1360 // Sum up fluxes
... ... @@ -1352,36 +1365,57 @@ double GSkymap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
1352 1365 } // endif: we had a HealPix map
1353 1366  
1354 1367 // ... otherwise perform flux computation for WCS map. The pixel is
1355   - // divided into 4 quadrants and the flux is computed in each quandrant
1356   - // by multiplying the intensity in the centre of the quadrant with the
1357   - // solid angle of the quadrant.
  1368 + // divided into 8 wedges and the flux is integrated in each wedge by
  1369 + // computing the average of the corner intensities multiplied by the
  1370 + // solid angle of the wedge. This leads to a precision of the order
  1371 + // of <1% in the flux computation.
1358 1372 else {
1359 1373  
1360   - // Compute flux
1361   - flux = this->operator()(pix2dir(GSkyPixel(pixel.x()-0.25,
1362   - pixel.y()-0.25)), map) *
1363   - solidangle(pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5)),
1364   - pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5)),
1365   - pix2dir(GSkyPixel(pixel.x(), pixel.y())),
1366   - pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()))) +
1367   - this->operator()(pix2dir(GSkyPixel(pixel.x()+0.25,
1368   - pixel.y()-0.25)), map) *
1369   - solidangle(pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5)),
1370   - pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5)),
1371   - pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y())),
1372   - pix2dir(GSkyPixel(pixel.x(), pixel.y()))) +
1373   - this->operator()(pix2dir(GSkyPixel(pixel.x()+0.25,
1374   - pixel.y()+0.25)), map) *
1375   - solidangle(pix2dir(GSkyPixel(pixel.x(), pixel.y())),
1376   - pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y())),
1377   - pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5)),
1378   - pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5))) +
1379   - this->operator()(pix2dir(GSkyPixel(pixel.x()-0.25,
1380   - pixel.y()+0.25)), map) *
1381   - solidangle(pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y())),
1382   - pix2dir(GSkyPixel(pixel.x(), pixel.y())),
1383   - pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5)),
1384   - pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5)));
  1374 + // Get centre
  1375 + GSkyDir centre = pix2dir(pixel);
  1376 +
  1377 + // Get boundaries
  1378 + GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
  1379 + GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
  1380 + GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
  1381 + GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()));
  1382 + GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
  1383 + GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5));
  1384 + GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
  1385 + GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
  1386 +
  1387 + // Compute intensities
  1388 + double i0 = this->operator()(centre);
  1389 + double i1 = this->operator()(boundary1);
  1390 + double i2 = this->operator()(boundary2);
  1391 + double i3 = this->operator()(boundary3);
  1392 + double i4 = this->operator()(boundary4);
  1393 + double i5 = this->operator()(boundary5);
  1394 + double i6 = this->operator()(boundary6);
  1395 + double i7 = this->operator()(boundary7);
  1396 + double i8 = this->operator()(boundary8);
  1397 +
  1398 + // Compute fluxes
  1399 + double flux1 = gammalib::onethird * (i1 + i2 + i0) *
  1400 + solidangle(boundary1, boundary2, centre);
  1401 + double flux2 = gammalib::onethird * (i2 + i3 + i0) *
  1402 + solidangle(boundary2, boundary3, centre);
  1403 + double flux3 = gammalib::onethird * (i3 + i4 + i0) *
  1404 + solidangle(boundary3, boundary4, centre);
  1405 + double flux4 = gammalib::onethird * (i4 + i5 + i0) *
  1406 + solidangle(boundary4, boundary5, centre);
  1407 + double flux5 = gammalib::onethird * (i5 + i6 + i0) *
  1408 + solidangle(boundary5, boundary6, centre);
  1409 + double flux6 = gammalib::onethird * (i6 + i7 + i0) *
  1410 + solidangle(boundary6, boundary7, centre);
  1411 + double flux7 = gammalib::onethird * (i7 + i8 + i0) *
  1412 + solidangle(boundary7, boundary8, centre);
  1413 + double flux8 = gammalib::onethird * (i8 + i1 + i0) *
  1414 + solidangle(boundary8, boundary1, centre);
  1415 +
  1416 + // Sum up fluxes
  1417 + flux = (flux1 + flux2 + flux3 + flux4 +
  1418 + flux5 + flux6 + flux7 + flux8);
1385 1419  
1386 1420 } // endelse: we had a WCS map
1387 1421  
... ... @@ -2602,7 +2636,7 @@ double GSkymap::solidangle(const GSkyDir&amp; dir1, const GSkyDir&amp; dir2,
2602 2636 std::tan(0.5*(s-a12)) *
2603 2637 std::tan(0.5*(s-a23)) *
2604 2638 std::tan(0.5*(s-a13))));
2605   -
  2639 +
2606 2640 // Return solid angle
2607 2641 return solidangle;
2608 2642 }
... ...