Sergio Hernández Cadena / gammalib

Commit 07220862fbe41cc3e20bbea8e946a7cf3631231c

Authored by Jürgen Knödlseder
1 parent f42b94b7

Add GSkyMap::region_circle() method (#3145)

Inline Side-by-side
Showing 5 changed files with 319 additions and 55 deletions
1 1 New Features and Important Changes in GammaLib 1.7.0
2 2  
3   -29 June 2020
  3 +30 June 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -44,6 +44,7 @@ The following methods have been added:
44 44 - GSkyDir::rotate() - Rotate sky direction (radians variant)
45 45 - GSkyDir::GSkyDir(GVector&) - Sky direction vector constructor
46 46 - GSkyMap::extract() - Extract section of sky map
  47 +- GSkyMap::region_circle() - Get bounding region circle of sky map
47 48 - GSkyRegionCircle::operator==() - Equality operator
48 49 - GSkyRegionCircle::operator!=() - Non-equality operator
49 50 - GModelPar::has_factor_min() - Test if parameter has factor minimum
... ... @@ -194,6 +195,9 @@ for the coordinate system that is specified by the coordsys argument. By default
194 195 this argument is set to "CEL", hence by default the methods return the position
195 196 angle with respect to celestial coordinates (#3199).
196 197  
  198 +Added GSkyMap::region_circle() method to compute a bounding region circle that
  199 +encloses all sky map pixels (#3145).
  200 +
197 201 The GSkyMap operators used the sky map value access operator that bi-linearly
198 202 interpolates sky map values even for maps with identical definition. This could
199 203 lead to rounding errors. For maps with identical definition the operators now
... ...
include/GSkyMap.hpp
  View file @0722086
1 1 /***************************************************************************
2 2 * GSkyMap.hpp - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -177,6 +177,7 @@ public:
177 177 const int& starty, const int& stopy) const;
178 178 GSkyMap extract(const GSkyRegions& inclusions) const;
179 179 void stack_maps(void);
  180 + GSkyRegionCircle region_circle(void) const;
180 181 void load(const GFilename& filename);
181 182 void save(const GFilename& filename,
182 183 const bool& clobber = false) const;
... ...
pyext/GSkyMap.i
  View file @0722086
1 1 /***************************************************************************
2 2 * GSkyMap.i - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -182,6 +182,7 @@ public:
182 182 const int& starty, const int& stopy) const;
183 183 GSkyMap extract(const GSkyRegions& inclusions) const;
184 184 void stack_maps(void);
  185 + GSkyRegionCircle region_circle(void) const;
185 186 void load(const GFilename& filename);
186 187 void save(const GFilename& filename,
187 188 const bool& clobber = false) const;
... ...
src/sky/GSkyMap.cpp
  View file @0722086
1 1 /***************************************************************************
2 2 * GSkyMap.cpp - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -2267,6 +2267,122 @@ void GSkyMap::stack_maps(void)
2267 2267  
2268 2268  
2269 2269 /***********************************************************************//**
  2270 + * @brief Return sky region circle that encloses the sky map
  2271 + *
  2272 + * @return Enclosing sky region circle.
  2273 + *
  2274 + * Returns a sky region circle that encloses the sky map.
  2275 + *
  2276 + * For a sky map in HealPix projection the method returns a region circle
  2277 + * that encloses the full sky, with a centre at Right Ascension and
  2278 + * Declination of 0 and a radius of 180 degrees.
  2279 + *
  2280 + * For a sky map in World Coordinate projection the method returns a region
  2281 + * circle that is centred on the central pixel of the sky map (i.e. the
  2282 + * pixel with indices \f$(n_x/2,n_y/2)\f$, where \f$n_x\f$ and \f$n_y\f$ are
  2283 + * the number of pixels in the X- and Y-drection) and a radius that is
  2284 + * determined from the maximum radial distance of all pixels with respect
  2285 + * to the central pixel.
  2286 + ***************************************************************************/
  2287 +GSkyRegionCircle GSkyMap::region_circle(void) const
  2288 +{
  2289 + // Initialise sky region circle
  2290 + GSkyRegionCircle region;
  2291 +
  2292 + // If the sky map is in HealPix projection then set radius to 180 deg
  2293 + if (projection()->code() == "HPX") {
  2294 + region = GSkyRegionCircle(0.0, 0.0, 180.0);
  2295 + }
  2296 +
  2297 + // ... otherwise sky map is in World Coordinate Projection
  2298 + else {
  2299 +
  2300 + // Initialise maximum radius
  2301 + double max_radius = 0.0;
  2302 + int max_index = -1;
  2303 +
  2304 + // Determine map centre
  2305 + GSkyPixel pixel(nx()/2.0, ny()/2.0);
  2306 + GSkyDir centre = pix2dir(pixel);
  2307 +
  2308 + // Determine maximum pixel distance with respect to map centre
  2309 + for (int i = 0; i < npix(); ++i) {
  2310 + double radius = inx2dir(i).dist_deg(centre);
  2311 + if (radius > max_radius) {
  2312 + max_radius = radius;
  2313 + max_index = i;
  2314 + }
  2315 + }
  2316 +
  2317 + // Consider edges of pixel with maximum distance
  2318 + if (max_index != -1) {
  2319 +
  2320 + // Set upper pixel edge just beyond the upper boundary since the
  2321 + // upper boundary is excluded
  2322 + const double up = 0.4999999;
  2323 +
  2324 + // Get pixel with maximum distance
  2325 + GSkyPixel pixel = inx2pix(max_index);
  2326 +
  2327 + // Get boundaries of pixel
  2328 + GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
  2329 + GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
  2330 + GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()-0.5));
  2331 + GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()));
  2332 + GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()+up));
  2333 + GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+up));
  2334 + GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+up));
  2335 + GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
  2336 +
  2337 + // Get radii
  2338 + double radius1 = boundary1.dist_deg(centre);
  2339 + double radius2 = boundary2.dist_deg(centre);
  2340 + double radius3 = boundary3.dist_deg(centre);
  2341 + double radius4 = boundary4.dist_deg(centre);
  2342 + double radius5 = boundary5.dist_deg(centre);
  2343 + double radius6 = boundary6.dist_deg(centre);
  2344 + double radius7 = boundary7.dist_deg(centre);
  2345 + double radius8 = boundary8.dist_deg(centre);
  2346 +
  2347 + // Determine maximum radius
  2348 + if (radius1 > max_radius) {
  2349 + max_radius = radius1;
  2350 + }
  2351 + if (radius2 > max_radius) {
  2352 + max_radius = radius2;
  2353 + }
  2354 + if (radius3 > max_radius) {
  2355 + max_radius = radius3;
  2356 + }
  2357 + if (radius4 > max_radius) {
  2358 + max_radius = radius4;
  2359 + }
  2360 + if (radius5 > max_radius) {
  2361 + max_radius = radius5;
  2362 + }
  2363 + if (radius6 > max_radius) {
  2364 + max_radius = radius6;
  2365 + }
  2366 + if (radius7 > max_radius) {
  2367 + max_radius = radius7;
  2368 + }
  2369 + if (radius8 > max_radius) {
  2370 + max_radius = radius8;
  2371 + }
  2372 +
  2373 + } // endif: considered edges of pixel
  2374 +
  2375 + // Set sky region
  2376 + region = GSkyRegionCircle(centre, max_radius);
  2377 +
  2378 + } // endelse: sky map was in World Coordinate Projection
  2379 +
  2380 + // Return region
  2381 + return region;
  2382 +}
  2383 +
  2384 +
  2385 +/***********************************************************************//**
2270 2386 * @brief Load skymap from FITS file.
2271 2387 *
2272 2388 * @param[in] filename FITS file name..
... ...
test/test_GSky.cpp
  View file @0722086
... ... @@ -578,60 +578,189 @@ void TestGSky::test_GSkyMap_healpix_construct(void)
578 578 {
579 579  
580 580 // Test void constructor
581   - test_try("Test void constructor");
582   - try {
583   - GSkyMap map;
584   - test_try_success();
585   - }
586   - catch (std::exception &e) {
587   - test_try_failure(e);
588   - }
  581 + GSkyMap empty_map;
  582 + test_assert(empty_map.is_empty(), "Check for empty sky map");
  583 + test_value(empty_map.npix(), 0, "Check for no sky map pixels");
  584 + test_value(empty_map.nx(), 0, "Check for no X-direction map pixels");
  585 + test_value(empty_map.ny(), 0, "Check for no Y-direction map pixels");
  586 + test_value(empty_map.nmaps(), 0, "Check for no sky maps");
589 587  
590   - // Test correct Healpix constructors
591   - test_try("Test correct Healpix constructors");
592   - try {
593   - GSkyMap ring1("GAL", 1, "RING", 1);
594   - GSkyMap ring2("GAL", 2, "RING", 1);
595   - GSkyMap ring4("GAL", 4, "RING", 1);
596   - GSkyMap ring8("GAL", 8, "RING", 1);
597   - GSkyMap ring16("GAL", 16, "RING", 1);
598   - GSkyMap ring32("GAL", 32, "RING", 1);
599   - GSkyMap ring64("GAL", 64, "RING", 1);
600   - GSkyMap ring128("GAL", 128, "RING", 1);
601   - GSkyMap ring256("GAL", 256, "RING", 1);
602   - GSkyMap ring512("GAL", 512, "RING", 1);
603   - GSkyMap nest1("GAL", 1, "NEST", 1);
604   - GSkyMap nest2("GAL", 2, "NEST", 1);
605   - GSkyMap nest4("GAL", 4, "NEST", 1);
606   - GSkyMap nest8("GAL", 8, "NEST", 1);
607   - GSkyMap nest16("GAL", 16, "NEST", 1);
608   - GSkyMap nest32("GAL", 32, "NEST", 1);
609   - GSkyMap nest64("GAL", 64, "NEST", 1);
610   - GSkyMap nest128("GAL", 128, "NEST", 1);
611   - GSkyMap nest256("GAL", 256, "NEST", 1);
612   - GSkyMap nest512("GAL", 512, "NEST", 1);
613   - GSkyMap map1("CEL", 1, "RING", 1);
614   - GSkyMap map2("CEL", 1, "RING", 2);
615   - GSkyMap map3("EQU", 1, "RING", 1);
616   - GSkyMap map4("EQU", 1, "NESTED", 1);
617   -
618   - test_try_success();
619   - }
620   - catch (std::exception &e) {
621   - test_try_failure(e);
622   - }
  588 + // Test Healpix constructors
  589 + GSkyMap ring1("GAL", 1, "RING", 1);
  590 + test_assert(!ring1.is_empty(), "Check for non-empty sky map");
  591 + test_value(ring1.npix(), 12, "Check for 12 sky map pixels");
  592 + test_value(ring1.nx(), 0, "Check for no X-direction map pixels");
  593 + test_value(ring1.ny(), 0, "Check for no Y-direction map pixels");
  594 + test_value(ring1.nmaps(), 1, "Check for one sky map");
  595 + //
  596 + GSkyMap ring2("GAL", 2, "RING", 1);
  597 + test_assert(!ring2.is_empty(), "Check for non-empty sky map");
  598 + test_value(ring2.npix(), 48, "Check for 48 sky map pixels");
  599 + test_value(ring2.nx(), 0, "Check for no X-direction map pixels");
  600 + test_value(ring2.ny(), 0, "Check for no Y-direction map pixels");
  601 + test_value(ring2.nmaps(), 1, "Check for one sky map");
  602 + //
  603 + GSkyMap ring4("GAL", 4, "RING", 1);
  604 + test_assert(!ring4.is_empty(), "Check for non-empty sky map");
  605 + test_value(ring4.npix(), 192, "Check for 192 sky map pixels");
  606 + test_value(ring4.nx(), 0, "Check for no X-direction map pixels");
  607 + test_value(ring4.ny(), 0, "Check for no Y-direction map pixels");
  608 + test_value(ring4.nmaps(), 1, "Check for one sky map");
  609 + //
  610 + GSkyMap ring8("GAL", 8, "RING", 1);
  611 + test_assert(!ring8.is_empty(), "Check for non-empty sky map");
  612 + test_value(ring8.npix(), 768, "Check for 768 sky map pixels");
  613 + test_value(ring8.nx(), 0, "Check for no X-direction map pixels");
  614 + test_value(ring8.ny(), 0, "Check for no Y-direction map pixels");
  615 + test_value(ring8.nmaps(), 1, "Check for one sky map");
  616 + //
  617 + GSkyMap ring16("GAL", 16, "RING", 1);
  618 + test_assert(!ring16.is_empty(), "Check for non-empty sky map");
  619 + test_value(ring16.npix(), 3072, "Check for 3072 sky map pixels");
  620 + test_value(ring16.nx(), 0, "Check for no X-direction map pixels");
  621 + test_value(ring16.ny(), 0, "Check for no Y-direction map pixels");
  622 + test_value(ring16.nmaps(), 1, "Check for one sky map");
  623 + //
  624 + GSkyMap ring32("GAL", 32, "RING", 1);
  625 + test_assert(!ring32.is_empty(), "Check for non-empty sky map");
  626 + test_value(ring32.npix(), 12288, "Check for 12288 sky map pixels");
  627 + test_value(ring32.nx(), 0, "Check for no X-direction map pixels");
  628 + test_value(ring32.ny(), 0, "Check for no Y-direction map pixels");
  629 + test_value(ring32.nmaps(), 1, "Check for one sky map");
  630 + //
  631 + GSkyMap ring64("GAL", 64, "RING", 1);
  632 + test_assert(!ring64.is_empty(), "Check for non-empty sky map");
  633 + test_value(ring64.npix(), 49152, "Check for 49152 sky map pixels");
  634 + test_value(ring64.nx(), 0, "Check for no X-direction map pixels");
  635 + test_value(ring64.ny(), 0, "Check for no Y-direction map pixels");
  636 + test_value(ring64.nmaps(), 1, "Check for one sky map");
  637 + //
  638 + GSkyMap ring128("GAL", 128, "RING", 1);
  639 + test_assert(!ring128.is_empty(), "Check for non-empty sky map");
  640 + test_value(ring128.npix(), 196608, "Check for 196608 sky map pixels");
  641 + test_value(ring128.nx(), 0, "Check for no X-direction map pixels");
  642 + test_value(ring128.ny(), 0, "Check for no Y-direction map pixels");
  643 + test_value(ring128.nmaps(), 1, "Check for one sky map");
  644 + //
  645 + GSkyMap ring256("GAL", 256, "RING", 1);
  646 + test_assert(!ring256.is_empty(), "Check for non-empty sky map");
  647 + test_value(ring256.npix(), 786432, "Check 786432 XX sky map pixels");
  648 + test_value(ring256.nx(), 0, "Check for no X-direction map pixels");
  649 + test_value(ring256.ny(), 0, "Check for no Y-direction map pixels");
  650 + test_value(ring256.nmaps(), 1, "Check for one sky map");
  651 + //
  652 + GSkyMap ring512("GAL", 512, "RING", 1);
  653 + test_assert(!ring512.is_empty(), "Check for non-empty sky map");
  654 + test_value(ring512.npix(), 3145728, "Check for 3145728 sky map pixels");
  655 + test_value(ring512.nx(), 0, "Check for no X-direction map pixels");
  656 + test_value(ring512.ny(), 0, "Check for no Y-direction map pixels");
  657 + test_value(ring512.nmaps(), 1, "Check for one sky map");
  658 + //
  659 + GSkyMap nest1("GAL", 1, "NEST", 1);
  660 + test_assert(!nest1.is_empty(), "Check for non-empty sky map");
  661 + test_value(nest1.npix(), 12, "Check for 12 sky map pixels");
  662 + test_value(nest1.nx(), 0, "Check for no X-direction map pixels");
  663 + test_value(nest1.ny(), 0, "Check for no Y-direction map pixels");
  664 + test_value(nest1.nmaps(), 1, "Check for one sky map");
  665 + //
  666 + GSkyMap nest2("GAL", 2, "NEST", 1);
  667 + test_assert(!nest2.is_empty(), "Check for non-empty sky map");
  668 + test_value(nest2.npix(), 48, "Check for 48 sky map pixels");
  669 + test_value(nest2.nx(), 0, "Check for no X-direction map pixels");
  670 + test_value(nest2.ny(), 0, "Check for no Y-direction map pixels");
  671 + test_value(nest2.nmaps(), 1, "Check for one sky map");
  672 + //
  673 + GSkyMap nest4("GAL", 4, "NEST", 1);
  674 + test_assert(!nest4.is_empty(), "Check for non-empty sky map");
  675 + test_value(nest4.npix(), 192, "Check for 192 sky map pixels");
  676 + test_value(nest4.nx(), 0, "Check for no X-direction map pixels");
  677 + test_value(nest4.ny(), 0, "Check for no Y-direction map pixels");
  678 + test_value(nest4.nmaps(), 1, "Check for one sky map");
  679 + //
  680 + GSkyMap nest8("GAL", 8, "NEST", 1);
  681 + test_assert(!nest8.is_empty(), "Check for non-empty sky map");
  682 + test_value(nest8.npix(), 768, "Check for 768 sky map pixels");
  683 + test_value(nest8.nx(), 0, "Check for no X-direction map pixels");
  684 + test_value(nest8.ny(), 0, "Check for no Y-direction map pixels");
  685 + test_value(nest8.nmaps(), 1, "Check for one sky map");
  686 + //
  687 + GSkyMap nest16("GAL", 16, "NEST", 1);
  688 + test_assert(!nest16.is_empty(), "Check for non-empty sky map");
  689 + test_value(nest16.npix(), 3072, "Check for 3072 sky map pixels");
  690 + test_value(nest16.nx(), 0, "Check for no X-direction map pixels");
  691 + test_value(nest16.ny(), 0, "Check for no Y-direction map pixels");
  692 + test_value(nest16.nmaps(), 1, "Check for one sky map");
  693 + //
  694 + GSkyMap nest32("GAL", 32, "NEST", 1);
  695 + test_assert(!nest32.is_empty(), "Check for non-empty sky map");
  696 + test_value(nest32.npix(), 12288, "Check for 12288 sky map pixels");
  697 + test_value(nest32.nx(), 0, "Check for no X-direction map pixels");
  698 + test_value(nest32.ny(), 0, "Check for no Y-direction map pixels");
  699 + test_value(nest32.nmaps(), 1, "Check for one sky map");
  700 + //
  701 + GSkyMap nest64("GAL", 64, "NEST", 1);
  702 + test_assert(!nest64.is_empty(), "Check for non-empty sky map");
  703 + test_value(nest64.npix(), 49152, "Check for 49152 sky map pixels");
  704 + test_value(nest64.nx(), 0, "Check for no X-direction map pixels");
  705 + test_value(nest64.ny(), 0, "Check for no Y-direction map pixels");
  706 + test_value(nest64.nmaps(), 1, "Check for one sky map");
  707 + //
  708 + GSkyMap nest128("GAL", 128, "NEST", 1);
  709 + test_assert(!nest128.is_empty(), "Check for non-empty sky map");
  710 + test_value(nest128.npix(), 196608, "Check for 196608 sky map pixels");
  711 + test_value(nest128.nx(), 0, "Check for no X-direction map pixels");
  712 + test_value(nest128.ny(), 0, "Check for no Y-direction map pixels");
  713 + test_value(nest128.nmaps(), 1, "Check for one sky map");
  714 + //
  715 + GSkyMap nest256("GAL", 256, "NEST", 1);
  716 + test_assert(!nest256.is_empty(), "Check for non-empty sky map");
  717 + test_value(nest256.npix(), 786432, "Check for 786432 sky map pixels");
  718 + test_value(nest256.nx(), 0, "Check for no X-direction map pixels");
  719 + test_value(nest256.ny(), 0, "Check for no Y-direction map pixels");
  720 + test_value(nest256.nmaps(), 1, "Check for one sky map");
  721 + //
  722 + GSkyMap nest512("GAL", 512, "NEST", 1);
  723 + test_assert(!nest512.is_empty(), "Check for non-empty sky map");
  724 + test_value(nest512.npix(), 3145728, "Check for 3145728 sky map pixels");
  725 + test_value(nest512.nx(), 0, "Check for no X-direction map pixels");
  726 + test_value(nest512.ny(), 0, "Check for no Y-direction map pixels");
  727 + test_value(nest512.nmaps(), 1, "Check for one sky map");
  728 + //
  729 + GSkyMap map1("CEL", 1, "RING", 1);
  730 + test_assert(!map1.is_empty(), "Check for non-empty sky map");
  731 + test_value(map1.npix(), 12, "Check for 12 sky map pixels");
  732 + test_value(map1.nx(), 0, "Check for no X-direction map pixels");
  733 + test_value(map1.ny(), 0, "Check for no Y-direction map pixels");
  734 + test_value(map1.nmaps(), 1, "Check for one sky map");
  735 + //
  736 + GSkyMap map2("CEL", 1, "RING", 2);
  737 + test_assert(!map2.is_empty(), "Check for non-empty sky map");
  738 + test_value(map2.npix(), 12, "Check for 12 sky map pixels");
  739 + test_value(map2.nx(), 0, "Check for no X-direction map pixels");
  740 + test_value(map2.ny(), 0, "Check for no Y-direction map pixels");
  741 + test_value(map2.nmaps(), 2, "Check for one sky map");
  742 + //
  743 + GSkyMap map3("EQU", 1, "RING", 1);
  744 + test_assert(!map3.is_empty(), "Check for non-empty sky map");
  745 + test_value(map3.npix(), 12, "Check for 12 sky map pixels");
  746 + test_value(map3.nx(), 0, "Check for no X-direction map pixels");
  747 + test_value(map3.ny(), 0, "Check for no Y-direction map pixels");
  748 + test_value(map3.nmaps(), 1, "Check for one sky map");
  749 + //
  750 + GSkyMap map4("EQU", 1, "NESTED", 2);
  751 + test_assert(!map4.is_empty(), "Check for non-empty sky map");
  752 + test_value(map4.npix(), 12, "Check for 12 sky map pixels");
  753 + test_value(map4.nx(), 0, "Check for no X-direction map pixels");
  754 + test_value(map4.ny(), 0, "Check for no Y-direction map pixels");
  755 + test_value(map4.nmaps(), 2, "Check for one sky map");
623 756  
624 757 // Test Healpix copy constructor
625   - test_try("Test Healpix copy constructor");
626   - try {
627   - GSkyMap ring1("GAL", 1, "RING", 1);
628   - GSkyMap ring2 = ring1;
629   -
630   - test_try_success();
631   - }
632   - catch (std::exception &e) {
633   - test_try_failure(e);
634   - }
  758 + GSkyMap ring1cpy = ring1;
  759 + test_assert(!ring1cpy.is_empty(), "Check for non-empty sky map");
  760 + test_value(ring1cpy.npix(), 12, "Check for 12 sky map pixels");
  761 + test_value(ring1cpy.nx(), 0, "Check for no X-direction map pixels");
  762 + test_value(ring1cpy.ny(), 0, "Check for no Y-direction map pixels");
  763 + test_value(ring1cpy.nmaps(), 1, "Check for one sky map");
635 764  
636 765 // Test invalid coordsys in constructor
637 766 test_try("Test invalid coordsys in constructor");
... ... @@ -686,6 +815,12 @@ void TestGSky::test_GSkyMap_healpix_construct(void)
686 815 test_try_failure(e);
687 816 }
688 817  
  818 + // Test region circle
  819 + GSkyRegionCircle region = ring1.region_circle();
  820 + test_value(region.centre().ra_deg(), 0.0, "Check R.A. for region_circle()");
  821 + test_value(region.centre().dec_deg(), 0.0, "Check Dec. for region_circle()");
  822 + test_value(region.radius(), 180.0, "Check radius for region_circle()");
  823 +
689 824 // Exit test
690 825 return;
691 826  
... ... @@ -813,6 +948,13 @@ void TestGSky::test_GSkyMap_wcs_construct(void)
813 948 test_try_failure(e);
814 949 }
815 950  
  951 + // Test region circle
  952 + GSkyMap map2("CAR", "CEL", 138.817, 37.293, 0.521, 0.931, 100, 100);
  953 + GSkyRegionCircle region = map2.region_circle();
  954 + test_value(region.centre().ra_deg(), 139.146485, "Check R.A. for region_circle()");
  955 + test_value(region.centre().dec_deg(), 37.758046, "Check Dec. for region_circle()");
  956 + test_value(region.radius(), 52.370103, "Check radius for region_circle()");
  957 +
816 958 // Exit test
817 959 return;
818 960 }
... ...