1!----------------------------------------------------------------------------- 2! 3! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, 4! Benjamin D. Wandelt, Anthony J. Banday, 5! Matthias Bartelmann, Hans K. Eriksen, 6! Frode K. Hansen, Martin Reinecke 7! 8! 9! This file is part of HEALPix. 10! 11! HEALPix is free software; you can redistribute it and/or modify 12! it under the terms of the GNU General Public License as published by 13! the Free Software Foundation; either version 2 of the License, or 14! (at your option) any later version. 15! 16! HEALPix is distributed in the hope that it will be useful, 17! but WITHOUT ANY WARRANTY; without even the implied warranty of 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19! GNU General Public License for more details. 20! 21! You should have received a copy of the GNU General Public License 22! along with HEALPix; if not, write to the Free Software 23! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 24! 25! For more information about HEALPix see http://healpix.sourceforge.net 26! 27!----------------------------------------------------------------------------- 28module pix_tools 29!================================================================== 30! Various subroutines for converting angles and unit vectors 31! to pixel numbers in both the NESTED and RING schemes, 32! as well as looking for neighbours and making pixel queries 33! 34! last update : Oct 4, 2002, EH, computation of pixel vertex in 35! pix2vec_ring and pix2vec_nest 36! 2004-05-28, EH: edition of next_in_line_nest for nside=1 37! addition of optional 'mask' in remove_dipole 38! 2004-06-01 correction of bug at low phi in in_ring 39! 2004-08-25 added multi-D convert_inplace_* and convert_nest2ring, convert_ring2nest 40! 2004-10-07 added template_pixel_*, same_shape_pixels_* 41! 2004-12-15 generic forms for convert_inplace (I,R,D with 1 or N dimensions), 42! remove_dipole (R,D) 43! 2005-01-25 ensure backward compatibility of remove_dipole 44! 2005-08-03 edit same_shape_pixels_nest to comply with GFORTRAN 45! 2006-06-28 edit to neighbours_nest to correct error in Nside=1 case 46! 2008-02-03: addition of weights in remove dipole 47! 2008-02-04: check ordering in [1,2] in remove dipole ; check nest in [0,1] in query_* 48! 2008-03-28: edit ring_num and query_strip to improve inclusive query_strip 49! bug fix on query_disc 50! 2009-03-04: starts supporting 64-bit npix 51! largest supported Nside is now 2^28 instead of 2^13 52! Note: the theoretical limit Nside=2^29 is not implemented to avoid 53! dealing with too many 64 bit integer variables that slow down execution 54! 2011-08-25 : 2011-08-*: improves accuracy of pixel routines close to Poles 55! 2011-10-13: improved query_triangle, introduced process_intervals 56! 2012-07-17: Parallel OpenMP implementation of medfiltmap 57! 2012-08-27: correction of a bug affecting neighbours_nest and next_in_line_nest at Nside>8192 58! 2013-04-02: bug correction in query_disc in inclusive mode 59! 2013-05-07: G95-compatible 60! 2015-09-02: added nest2uniq, uniq2nest 61! 2018-05-22: added nside2npweights 62! 2018-10-22: added apply_mask 63!================================================================== 64 ! subroutine query_strip Done (To be Tested) depends on in_ring 65 ! subroutine query_polygon Done (To be Tested) depends on isort 66 ! subroutine query_triangle Done (To be Tested) 67 ! subroutine getdisc_ring Done depends on query_disc 68 ! subroutine query_disc Done (To be Tested) 69 ! function ring_num OK 70 ! function ring2z Done 71 ! subroutine in_ring Done (To be Tested) depends in next_in_line_nest 72 ! subroutine intrs_intrv OK 73 ! subroutine process_intervals ok 74 ! function fudge_query_radius ok 75 !!!!!! subroutine correct_ring_phi OK 76 ! 77 ! subroutine pix2ang_ring Done 78 ! subroutine pix2vec_ring Done 79 ! subroutine ang2pix_ring Done 80 ! subroutine vec2pix_ring Done 81 ! subroutine ang2pix_nest Done 82 ! subroutine vec2pix_nest Done 83 ! subroutine nest2ring Done 84 ! subroutine ring2nest Done 85 ! function npix2nside Done 86 ! 87 ! subroutine warning_oldbounds OK 88 ! subroutine remove_dipole_real Done depends on query_strip 89 ! subroutine remove_dipole_double Done 90 ! subroutine remove_dipole_real_old Done 91 ! subroutine remove_dipole_double_old Done 92 ! subroutine remove_dipole_real_v12 Done 93 ! subroutine remove_dipole_double_v12 Done 94 ! subroutine add_dipole_real Done 95 ! subroutine add_dipole_double Done 96 ! subroutine medfiltmap_S Done (To be tested) 97 ! subroutine medfiltmap_D Done (To be tested) 98 99 ! subroutine convert_nest2ring_int_1d Done 100 ! subroutine convert_nest2ring_int8_1d Done 101 ! subroutine convert_nest2ring_real_1d Done 102 ! subroutine convert_nest2ring_double_1d Done 103 ! subroutine convert_ring2nest_int_1d Done 104 ! subroutine convert_ring2nest_int8_1d Done 105 ! subroutine convert_ring2nest_real_1d Done 106 ! subroutine convert_ring2nest_double_1d Done 107 ! subroutine convert_nest2ring_int_nd Done 108 ! subroutine convert_nest2ring_int8_nd Done 109 ! subroutine convert_nest2ring_real_nd Done 110 ! subroutine convert_nest2ring_double_nd Done 111 ! subroutine convert_ring2nest_int_nd Done 112 ! subroutine convert_ring2nest_int8_nd Done 113 ! subroutine convert_ring2nest_real_nd Done 114 ! subroutine convert_ring2nest_double_nd Done 115 ! 116 ! subroutine convert_inplace_int_1d Done (to be tested) 117 ! subroutine convert_inplace_real_1d Done (to be tested) 118 ! subroutine convert_inplace_double_1d Done (to be tested) 119 ! subroutine convert_inplace_int_nd Done (to be tested) 120 ! subroutine convert_inplace_real_nd Done (to be tested) 121 ! subroutine convert_inplace_double_nd Done (to be tested) 122 ! 123 ! subroutine xy2pix_nest Done (to be tested) 124 ! subroutine pix2xy_nest Done (to be tested) 125 ! subroutine mk_pix2xy OK 126 ! subroutine mk_xy2pix OK 127 ! subroutine mk_xy2pix1 OK 128 ! subroutine neighbours_nest Done (to be tested) 129 ! subroutine next_in_line_nest Done depends on xy2pix_nest pix2xy_nest 130 ! subroutine ang2vec OK 131 ! subroutine vec2ang OK 132 ! function nside2npix Done 133 ! subroutine surface_triangle OK 134 ! subroutine angdist OK 135 ! subroutine vect_prod OK 136 ! function nside2ntemplates Done 137 ! subroutine template_pixel_ring Done (to be tested) 138 ! subroutine same_shape_pixels_ring Done (to be tested) 139 ! subroutine template_pixel_nest Done (to be tested) 140 ! subroutine same_shape_pixels_nest Done (to be tested) 141! 142! 2011-08-*: improves accuracy close to Poles 143! ang2pix_nest small 144! ang2pix_ring small 145! ang2vec OK 146! angdist Done 147! pix2ang_nest small, kshift 148! pix2ang_ring isqrt, small 149! pix2vec_nest small, [kshift TODO] 150! pix2vec_ring isqrt, small 151! cheap_isqrt Added 152! vec2ang Done 153! vec2pix_nest small 154! vec2pix_ring small 155! ring2nest isqrt 156 157 USE healpix_types 158 USE misc_utils 159 USE long_intrinsic, only: long_size 160 IMPLICIT none 161 162 INTEGER(KIND=i4b), private, PARAMETER :: ns_max4=8192 ! 2^13 163 INTEGER(KIND=i4b), private, PARAMETER :: ns_max8=268435456! 2^28 164#ifdef NO64BITS 165 INTEGER(KIND=i4b), private, PARAMETER :: ns_max=ns_max4 ! largest nside available 166#else 167 INTEGER(KIND=i4b), private, PARAMETER :: ns_max=ns_max8 ! largest nside available 168#endif 169 170 !initialise array x2pix, y2pix and pix2x, pix2y used in several routines 171! integer(KIND=i4b), private, save, dimension(128) :: x2pix=0,y2pix=0 172 integer(KIND=i4b), private, save, dimension(0:127) :: x2pix1=-1,y2pix1=-1 173 174 integer(KIND=i4b), private, save, dimension(0:1023) :: pix2x=-1, pix2y=-1 175 176 ! coordinate of the lowest corner of each face 177 INTEGER(KIND=I4B), private, save, dimension(0:11) :: jrll1 = (/ 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 /) ! in unit of nside 178 INTEGER(KIND=I4B), private, save, dimension(0:11) :: jpll1 = (/ 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 /) ! in unit of nside/2 179 180 ! obsolete 181 interface convert_inplace_real 182 module procedure convert_inplace_real_1d, convert_inplace_real_nd 183 end interface 184 interface convert_inplace_int 185 module procedure convert_inplace_int_1d, convert_inplace_int_nd 186 end interface 187 188 ! generic form 189 interface convert_inplace 190 module procedure convert_inplace_real_1d, convert_inplace_real_nd, & 191 & convert_inplace_int_1d, convert_inplace_int_nd, & 192 & convert_inplace_double_1d, convert_inplace_double_nd 193 end interface 194 195 interface remove_dipole 196 module procedure remove_dipole_real, remove_dipole_double, & 197 & remove_dipole_real_old, remove_dipole_double_old, & 198 & remove_dipole_real_v12, remove_dipole_double_v12 199 end interface 200 201 interface add_dipole 202 module procedure add_dipole_real, add_dipole_double 203 end interface 204 205 interface apply_mask 206 module procedure apply_mask_real, apply_mask_double 207 end interface 208 209 interface convert_nest2ring 210 module procedure convert_nest2ring_int_1d, & 211 & convert_nest2ring_real_1d, & 212 & convert_nest2ring_double_1d, & 213 & convert_nest2ring_int_nd, & 214 & convert_nest2ring_real_nd, & 215 & convert_nest2ring_double_nd 216#ifndef NO64BITS 217 module procedure convert_nest2ring_int8_1d, & 218 convert_nest2ring_int8_nd 219#endif 220 end interface 221 222 interface convert_ring2nest 223 module procedure convert_ring2nest_int_1d, & 224 & convert_ring2nest_real_1d, & 225 & convert_ring2nest_double_1d, & 226 & convert_ring2nest_int_nd, & 227 & convert_ring2nest_real_nd, & 228 & convert_ring2nest_double_nd 229#ifndef NO64BITS 230 module procedure convert_ring2nest_int8_1d, & 231 convert_ring2nest_int8_nd 232#endif 233 end interface 234 235 interface medfiltmap 236 module procedure medfiltmap_s, medfiltmap_d 237 end interface 238 239#ifdef NO64BITS 240 interface npix2nside 241 module procedure npix2nside 242 end interface 243 interface ang2pix_nest 244 module procedure ang2pix_nest 245 end interface 246 interface vec2pix_nest 247 module procedure vec2pix_nest 248 end interface 249 interface ang2pix_ring 250 module procedure ang2pix_ring 251 end interface 252 interface vec2pix_ring 253 module procedure vec2pix_ring 254 end interface 255 interface ring2nest 256 module procedure ring2nest 257 end interface 258 interface nest2ring 259 module procedure nest2ring 260 end interface 261 interface pix2ang_ring 262 module procedure pix2ang_ring 263 end interface 264 interface pix2vec_ring 265 module procedure pix2vec_ring 266 end interface 267 interface pix2ang_nest 268 module procedure pix2ang_nest 269 end interface 270 interface pix2vec_nest 271 module procedure pix2vec_nest 272 end interface 273 interface in_ring 274 module procedure in_ring 275 end interface 276 interface getdisc_ring 277 module procedure getdisc_ring 278 end interface 279 interface query_strip 280 module procedure query_strip 281 end interface 282 interface query_disc_old ! <<<<<<<<<<<<<<<<<<<<<< 283 module procedure query_disc_old 284 end interface 285 interface query_disc 286 module procedure query_disc 287 end interface 288 interface query_polygon 289 module procedure query_polygon 290 end interface 291 interface query_triangle 292 module procedure query_triangle 293 end interface 294 interface xy2pix_nest 295 module procedure xy2pix_nest 296 end interface 297 interface pix2xy_nest 298 module procedure pix2xy_nest 299 end interface 300 interface neighbours_nest 301 module procedure neighbours_nest 302 end interface 303 interface next_in_line_nest 304 module procedure next_in_line_nest 305 end interface 306 interface template_pixel_ring 307 module procedure template_pixel_ring 308 end interface 309 interface same_shape_pixels_ring 310 module procedure same_shape_pixels_ring 311 end interface 312 interface template_pixel_nest 313 module procedure template_pixel_nest 314 end interface 315 interface same_shape_pixels_nest 316 module procedure same_shape_pixels_nest 317 end interface 318 interface cheap_isqrt 319 module procedure cheap_isqrt_4 320 end interface 321 interface discedge2fulldisc 322 module procedure discedge2fulldisc 323 end interface 324 interface nest2uniq 325 module procedure nest2uniq 326 end interface 327 interface uniq2nest 328 module procedure uniq2nest 329 end interface 330#else 331 interface npix2nside 332 module procedure npix2nside, npix2nside_8 333 end interface 334 interface ang2pix_nest 335 module procedure ang2pix_nest, ang2pix_nest_8 336 end interface 337 interface vec2pix_nest 338 module procedure vec2pix_nest, vec2pix_nest_8 339 end interface 340 interface ang2pix_ring 341 module procedure ang2pix_ring, ang2pix_ring_8 342 end interface 343 interface vec2pix_ring 344 module procedure vec2pix_ring, vec2pix_ring_8 345 end interface 346 interface ring2nest 347 module procedure ring2nest, ring2nest_8 348 end interface 349 interface nest2ring 350 module procedure nest2ring, nest2ring_8 351 end interface 352 interface pix2ang_ring 353 module procedure pix2ang_ring, pix2ang_ring_8 354 end interface 355 interface pix2vec_ring 356 module procedure pix2vec_ring, pix2vec_ring_8 357 end interface 358 interface pix2ang_nest 359 module procedure pix2ang_nest, pix2ang_nest_8 360 end interface 361 interface pix2vec_nest 362 module procedure pix2vec_nest, pix2vec_nest_8 363 end interface 364 interface in_ring 365 module procedure in_ring, in_ring_8 366 end interface 367 interface getdisc_ring 368 module procedure getdisc_ring, getdisc_ring_8 369 end interface 370 interface query_strip 371 module procedure query_strip, query_strip_8 372 end interface 373 interface query_disc_old ! <<<<<<<<<<<<<<<<<<<<<< 374 module procedure query_disc_old, query_disc_old_8 375 end interface 376 interface query_disc 377 module procedure query_disc, query_disc_8 378 end interface 379 interface query_polygon 380 module procedure query_polygon, query_polygon_8 381 end interface 382 interface query_triangle 383 module procedure query_triangle, query_triangle_8 384 end interface 385 interface xy2pix_nest 386 module procedure xy2pix_nest, xy2pix_nest_8 387 end interface 388 interface pix2xy_nest 389 module procedure pix2xy_nest, pix2xy_nest_8 390 end interface 391 interface neighbours_nest 392 module procedure neighbours_nest, neighbours_nest_8 393 end interface 394 interface next_in_line_nest 395 module procedure next_in_line_nest, next_in_line_nest_8 396 end interface 397 interface template_pixel_ring 398 module procedure template_pixel_ring, template_pixel_ring_8 399 end interface 400 interface same_shape_pixels_ring 401 module procedure same_shape_pixels_ring, same_shape_pixels_ring_8 402 end interface 403 interface template_pixel_nest 404 module procedure template_pixel_nest, template_pixel_nest_8 405 end interface 406 interface same_shape_pixels_nest 407 module procedure same_shape_pixels_nest, same_shape_pixels_nest_8 408 end interface 409 interface cheap_isqrt 410 module procedure cheap_isqrt_4, cheap_isqrt_8 411 end interface 412 interface discedge2fulldisc 413 module procedure discedge2fulldisc, discedge2fulldisc_8 414 end interface 415 interface nest2uniq 416 module procedure nest2uniq, nest2uniq_8 417 end interface 418 interface uniq2nest 419 module procedure uniq2nest, uniq2nest_8 420 end interface 421#endif 422 423 private 424 425 public :: remove_dipole, add_dipole, & 426 & query_strip, & 427 & query_polygon, & 428 & query_triangle, & 429 & query_disc, & 430 & pix2ang_ring, pix2vec_ring, ang2pix_ring, vec2pix_ring, & 431 & pix2ang_nest, pix2vec_nest, ang2pix_nest, vec2pix_nest, & 432 & convert_nest2ring, convert_ring2nest, & 433 & convert_inplace, convert_inplace_real, convert_inplace_int, & 434 & nest2ring, ring2nest, xy2pix_nest, pix2xy_nest, & 435!!! & mk_pix2xy, mk_xy2pix, & 436 & neighbours_nest, & 437 & next_in_line_nest, & 438 & ang2vec, vec2ang, & 439 & npix2nside, nside2npix, & 440 & surface_triangle, angdist, vect_prod 441 442! public :: ang2pix_nest_old, vec2pix_nest_old, & 443! ang2pix_ring_old, vec2pix_ring_old, & 444! ring2nest_old, nest2ring_old, & 445! pix2ang_ring_old, pix2vec_ring_old, & 446! pix2ang_nest_old, pix2vec_nest_old, & 447! template_pixel_ring_old, same_shape_pixels_ring_old, & 448! template_pixel_nest_old, same_shape_pixels_nest_old 449 450 public :: query_disc_old 451 452 public :: nside2ntemplates, & 453 & template_pixel_ring, same_shape_pixels_ring, & 454 & template_pixel_nest, same_shape_pixels_nest 455 456 public :: medfiltmap 457 458 public :: intrs_intrv, in_ring, ring_num, ring2z ! arcane usage 459 public :: process_intervals, fudge_query_radius, discphirange_at_z 460 461 public :: getdisc_ring ! obsolete 462 463 public :: nest2uniq, uniq2nest 464 465 public :: nside2npweights 466 467 public :: apply_mask 468 469 470contains 471 472 473! !======================================================================= 474! function discphirange_at_z (vcenter, radius, z, phi0) result(dphi) 475! !======================================================================= 476! ! for a disc centered on vcenter and given radius, 477! ! and for location z=cos(theta) 478! ! returns dphi such that the disc has abs(phi-phi0) <= dphi 479! ! 480! ! solving disc equation on sphere: 481! ! sin(theta) sin(theta0) cos(phi-phi0) + cos(theta) cos(theta0) >= cos(radius) 482! ! 483! ! 2011-06-21: adapted from IDL routine of same name 484! !======================================================================= 485! real(DP), dimension(1:), intent(in) :: vcenter 486! real(DP), intent(in) :: radius, z 487! real(DP), intent(out) :: phi0 488! real(DP) :: dphi 489! ! 490! real(DP) :: cosang, cosdphi, norm, x0, y0, z0, st0, diff, st 491! real(DP), parameter :: zero=0.0_dp, one=1.0_dp 492 493! cosang = cos(radius) 494! norm = sqrt(sum(vcenter(1:3)**2)) 495! x0 = vcenter(1) / norm 496! y0 = vcenter(2) / norm 497! z0 = vcenter(3) / norm 498 499! phi0 = zero 500! if ((x0 /= zero) .or. (y0 /= zero)) phi0 = atan2(y0, x0) 501 502! st0 = x0*x0 + y0*y0 ! sin(theta0)^2 503! diff = cosang - z*z0 ! cos(rad) - cos(theta)cos(theta0) 504 505! dphi = -1000.0_dp ! out of disc 506! if (st0 == zero) then ! polar cap 507! if (diff <= zero) dphi = PI ! fully in cap 508! else 509! st = max(one - z*z, 1.0e-12_dp)! sin(theta)^2 510! cosdphi = diff/sqrt(st0*st) 511! if (cosdphi < -one) dphi = PI ! fully in disc 512! if (abs(cosdphi) <= one) dphi = acos(cosdphi) 513! endif 514 515! end function discphirange_at_z 516 !======================================================================= 517 function ring_num (nside, z, shift) result(ring_num_result) 518 !======================================================================= 519 ! ring = ring_num(nside, z [, shift=]) 520 ! returns the ring number in {1, 4*nside-1} 521 ! from the z coordinate 522 ! usually returns the ring closest to the z provided 523 ! if shift < 0, returns the ring immediatly north (of smaller index) of z 524 ! if shift > 0, returns the ring immediatly south (of smaller index) of z 525 ! 526 !======================================================================= 527 INTEGER(KIND=I4B) :: ring_num_result 528 REAL(KIND=DP), INTENT(IN) :: z 529 INTEGER(KIND=I4B), INTENT(IN) :: nside 530 integer(i4b), intent(in), optional :: shift 531 532 INTEGER(KIND=I4B) :: iring 533 real(DP) :: my_shift 534 !======================================================================= 535 536 537 my_shift = 0.0_dp 538 if (present(shift)) my_shift = shift * 0.5_dp 539 540 ! ----- equatorial regime --------- 541 iring = NINT( nside*(2.0_dp-1.500_dp*z) + my_shift ) 542 543 ! ----- north cap ------ 544 if (z > twothird) then 545 iring = NINT( nside* SQRT(3.0_dp*(1.0_dp-z)) + my_shift ) 546 if (iring == 0) iring = 1 547 endif 548 549 ! ----- south cap ----- 550 if (z < -twothird ) then 551 ! beware that we do a -shift in the south cap 552 iring = NINT( nside* SQRT(3.0_dp*(1.0_dp+z)) - my_shift ) 553 if (iring == 0) iring = 1 554 iring = 4*nside - iring 555 endif 556 557 ring_num_result = iring 558 559 return 560 end function ring_num 561 !======================================================================= 562 function ring2z (nside, ir) result(z) 563 !======================================================================= 564 ! returns the z coordinate of ring ir for Nside 565 ! 2009-03-25: accepts Nside > 8192 566 !======================================================================= 567 integer(kind=I4B), intent(in) :: ir 568 integer(kind=I4B), intent(in) :: nside 569 real(kind=DP) :: z 570 571 real(DP) :: fn, tmp 572 !======================================================================= 573 574 fn = real(nside,kind=DP) 575 576 if (ir < nside) then ! polar cap (north) 577 tmp = real(ir, DP) 578 z = 1.0_dp - (tmp * tmp) / (3.0_dp * fn * fn) 579 else if (ir < 3*nside) then ! tropical band 580 z = real( 2*nside-ir, kind=DP) * 2.0_dp / (3.0_dp * fn) 581 else ! polar cap (south) 582 tmp = real(4*nside - ir, DP) 583 z = -1.0_dp + (tmp * tmp) / (3.0_dp * fn * fn) 584 endif 585 586 return 587 end function ring2z 588 !======================================================================= 589 subroutine intrs_intrv( d1, d2, di, ni) 590 !======================================================================= 591 ! computes the intersection di 592 ! of 2 intervals d1 (= [a1,b1]) and d2 (= [a2,b2]) 593 ! on the periodic domain ( = [A,B], where A and B are arbitrary) 594 ! ni is the resulting number of intervals (0,1, or 2) 595 ! 596 ! if a1<b1 then d1 = {x | a1 <= x <= b1} 597 ! if a1>b1 then d1 = {x | a1 <= x <= B U A <= x <= b1} 598 !======================================================================= 599 real(kind=DP), dimension(1:), INTENT(IN) :: d1, d2 600 real(kind=DP), dimension(1:), INTENT(OUT) :: di 601 integer(kind=I4B), INTENT(OUT) :: ni 602 603 real(kind=DP), dimension(1:4) :: dk 604 integer(kind=I4B) :: ik 605 logical(kind=LGT) :: tr12, tr21, tr34, tr43, tr13, tr31, tr24, tr42, tr14, tr32 606 !======================================================================= 607 608 tr12 = (d1(1) < d1(2)) 609 tr21 = .NOT. tr12 610 tr34 = (d2(1) < d2(2)) 611 tr43 = .NOT. tr34 612 tr13 = (d1(1) < d2(1)) 613 tr31 = .NOT. tr13 614 tr24 = (d1(2) < d2(2)) 615 tr42 = .NOT. tr24 616 tr14 = (d1(1) < d2(2)) 617 tr32 = (d2(1) < d1(2)) 618 619 ik = 0 620 dk(1:4) = -1.0e9_dp 621 622 623 if ((tr31.AND.tr14) .OR. (tr43.AND.(tr31.OR.tr14))) then 624 ik = ik + 1 625 dk(ik) = d1(1) ! a1 626 endif 627 if ((tr13.AND.tr32) .OR. (tr21.AND.(tr13.OR.tr32))) then 628 ik = ik + 1 629 dk(ik) = d2(1) ! a2 630 endif 631 if ((tr32.AND.tr24) .OR. (tr43.AND.(tr32.OR.tr24))) then 632 ik = ik + 1 633 dk(ik) = d1(2) ! b1 634 endif 635 if ((tr14.AND.tr42) .OR. (tr21.AND.(tr14.OR.tr42))) then 636 ik = ik + 1 637 dk(ik) = d2(2) ! b2 638 endif 639 640 di(1:4) = 0.0_dp 641 select case (ik) 642 case (0) 643 ni = 0 644 case (2) 645 ni = 1 646 di(1:2) = (/ dk(1), dk(2) /) ! [a1,b1] or [a1,b2] or [a2,b1] or [a2,b2] 647 case (4) 648 ni = 2 649 di(1:4) = (/ dk(1), dk(4), dk(2), dk(3) /) ! [a1,b2] U [a2,b1] 650 case default 651 print*,"error in intrs_intrv", ik 652 print*,dk 653 print*,d1,d2 654 end select 655 656 return 657 end subroutine intrs_intrv 658 659 !======================================================================= 660 subroutine process_intervals(interval1, interval_list, interval_out, n_out) 661 !======================================================================= 662 ! intersection of 1 interval (defined by its 2 boundaries) 663 ! with an arbitrary list of intervals (defined as n*2 boundaries) 664 !======================================================================= 665 real(dp), dimension(1:), intent(in) :: interval1 666 real(dp), dimension(1:), intent(in) :: interval_list 667 real(dp), dimension(1:), intent(out):: interval_out 668 integer(i4b), intent(out) :: n_out 669 ! 670 integer(i4b) :: n_in, i_in, n_tmp_out 671 real(dp), dimension(1:4) :: w4 672 real(dp), dimension(1:30) :: work 673 !======================================================================= 674 675 n_out = 0 676 interval_out = 0.d0 677 n_in = size(interval_list)/2 678 work(1:2*n_in) = interval_list(1:2*n_in) ! copy input array to avoid overwritting it 679 680 if (n_in > 0) then 681 do i_in=0, n_in-1 682 ! intersection of 2 intervals -> 0,1,2 intervals 683 call intrs_intrv(interval1, work(1+2*i_in:2+2*i_in), w4, n_tmp_out) 684 if (n_tmp_out > 0) then ! n_tmp_out = 1 or 2 685 interval_out(2*n_out+1: 2*(n_out+n_tmp_out)) = w4(1:2*n_tmp_out) 686 n_out = n_out + n_tmp_out 687 endif 688 enddo 689 endif 690 691 return 692 end subroutine process_intervals 693 !======================================================================= 694 function fudge_query_radius(nside, radius_in, quadratic) result(radius_out) 695 !======================================================================= 696 ! radius_out = fudge_query_radius( nside, radius_in, QUADRATIC=) 697 ! 698 ! with 699 ! radius_out = radius_in + fudge (default) 700 ! or 701 ! radius_out = sqrt(radius_in^2 + fudge^2) (quadratic) 702 ! 703 ! if absent, radius_in = 0 704 ! where 705 ! fudge = factor(nside) * Pi / (4 *Nside) 706 ! with factor = sqrt( 5/9 + (8/Pi)^2 / 5 * (1-1/(2*Nside)) ) 707 ! 708 ! an upper bound of the actual largest radius 709 ! (determined by pixel located at z=2/3 and its North corner). 710 ! 711 ! 712 ! 2011-10-14, EH, v1.0 713 !======================================================================= 714 integer(i4b), intent(in) :: nside 715 real(dp), intent(in), optional :: radius_in 716 logical(lgt), intent(in), optional :: quadratic 717 real(dp) :: radius_out 718 logical(LGT) :: do_quad 719 real(dp) :: factor, fudge 720 !======================================================================= 721 radius_out = 0.0_dp 722 if (present(radius_in)) radius_out = radius_in 723 do_quad = .false. 724 if (present(quadratic)) do_quad = quadratic 725 726 factor = sqrt( 0.55555555555555555555_dp + 1.29691115062192347448165712908_dp*(1.0_dp-0.5_dp/nside) ) 727 fudge = factor * PI / (4.0_dp * nside) ! increase effective radius 728 729 if (do_quad) then 730 radius_out = sqrt(radius_out**2 + fudge**2) 731 else 732 radius_out = radius_out + fudge 733 endif 734 radius_out = min(radius_out, PI) 735 736 return 737 end function fudge_query_radius 738 !======================================================================= 739 subroutine discphirange_at_z(vector0, radius, z, nz, dphi, phi0) 740 !======================================================================= 741 ! for a disc centered on vcenter and given radius, 742 ! and for location z=cos(theta) 743 ! returns dphi such that the disc has abs(phi-phi0) <= dphi 744 ! 745 ! solving disc equation on sphere: 746 ! sin(theta) sin(theta0) cos(phi-phi0) + cos(theta) cos(theta0) >= cos(radius) 747 ! 748 !======================================================================= 749 real(dp), dimension(1:3), intent(in) :: vector0 750 real(dp), intent(in) :: radius 751 real(dp), dimension(1:), intent(in) :: z 752 integer(i4b), intent(in) :: nz 753 real(dp), dimension(1:), intent(out) :: dphi 754 real(dp), intent(out) :: phi0 755 ! 756 real(dp) :: cosang, cosphi0, x0, y0, z0, a, b, c, cosdphi, norm 757 integer(i4b) :: i 758 real(DP), parameter :: zero=0.0_dp, one=1.0_dp 759 760 cosang = cos(radius) 761 norm = sqrt(sum(vector0(1:3)**2)) 762 x0 = vector0(1) / norm 763 y0 = vector0(2) / norm 764 z0 = vector0(3) / norm 765 766 phi0 = zero 767 if ((x0 /= zero).or.(y0 /= zero)) phi0 = atan2(y0, x0) ! in ]-Pi, Pi] 768 !cosphi0 = cos(phi0) 769 a = x0*x0 + y0*y0 ! sin(theta0)^2 770 771 do i=1, nz 772 dphi(i) = -1000.0_dp ! default value for out of disc 773 774 b = cosang - z(i)*z0 ! cos(rad) - cos(theta)cos(theta0) 775 if (a == zero) then ! poles 776 if (b <= zero) dphi(i) = PI 777 else 778 c = max(one - z(i)*z(i), 1.0e-12_dp) ! sin(theta)^2 779 cosdphi = b / sqrt(a*c) 780 if (cosdphi < -one) dphi(i) = PI ! all the pixels at this elevation are in the disc 781 if (abs(cosdphi) <= one) dphi(i) = acos(cosdphi) ! in [0,Pi] 782 endif 783 784 enddo 785 786 return 787 end subroutine discphirange_at_z 788 !======================================================================= 789 subroutine pixels_on_edge(nside, irmin, irmax, phi0, dphi, ringphi, ngr) 790 !======================================================================= 791 !======================================================================= 792 integer(i4b), intent(in) :: nside, irmin, irmax 793 real(dp), intent(in) :: phi0 794 real(dp), dimension(1:), intent(in) :: dphi 795 integer(i4b), dimension(1:,1:), intent(out) :: ringphi 796 integer(i4b), intent(out) :: ngr 797 integer(i4b), parameter :: badvalue = -9999 798 integer(i4b) :: nrings, ir, k, thisring, npr 799 real(dp) , parameter :: zero = 0.0_dp 800 integer(i4b) :: iphi_low, iphi_hi, kshift 801 real(dp) :: shift 802 803 nrings = irmax - irmin + 1 804 805 ngr = 0 806 do thisring = irmin, irmax 807 ir = thisring - irmin + 1 808 call pixels_per_ring(nside, thisring, npr, kshift) 809 if (dphi(ir) >= PI) then ! full ring 810 ngr = ngr + 1 811 ringphi(1, ngr) = thisring 812 ringphi(2, ngr) = 0 813 ringphi(3, ngr) = npr-1 814 elseif (dphi(ir) >= zero) then ! partial ring 815 shift = kshift * 0.5_dp 816 iphi_low = ceiling (npr * (phi0 - dphi(ir)) / TWOPI - shift) 817 iphi_hi = floor (npr * (phi0 + dphi(ir)) / TWOPI - shift) 818 if (iphi_hi >= iphi_low) then ! pixel center in range 819 ngr = ngr + 1 820 ringphi(1, ngr) = thisring 821 ringphi(2, ngr) = modulo(iphi_low, npr) 822 ringphi(3, ngr) = modulo(iphi_hi, npr) 823 endif 824 endif 825 enddo 826 827 end subroutine pixels_on_edge 828 !======================================================================= 829 subroutine pixels_per_ring(nside, ring, npr, kshift, npnorth) 830 !======================================================================= 831 ! for a given Nside and ring index in [1,4*Nside-1], 832 ! returns the number of pixels in ring, their shift (0 or 1) in azimuth 833 ! and the number of pixels in current ring and above (=North) 834 ! 835 ! NB: 'rings' 0 and 4*Nside respectively are North and South Poles 836 !======================================================================= 837 integer(i4b), intent(in) :: nside, ring 838 integer(i4b), intent(out) :: npr, kshift 839 integer(i8b), intent(out), optional :: npnorth 840 integer(i8b) :: ncap, npix, ir 841 842 ! number of pixels in current ring 843 npr = min(nside, ring, 4*nside-ring) * 4 844 ! shift 845 kshift = mod(ring + 1, 2) ! 1 for even, 0 for odd 846 if (nside == 1) kshift = 1 - kshift ! except for Nside=1 847 if (npr < 4*nside) kshift = 1 ! 1 on polar cap 848 ! Number of pixels in current ring and above 849 if (present(npnorth)) then 850 if (ring <= nside) then ! in North cap 851 npnorth = ring*(ring+1_i8b)*2_i8b 852 elseif (ring <= 3*nside) then ! in Equatorial region 853 ncap = nside*(nside+1_i8b)*2_i8b 854 ir = ring-nside 855 npnorth = ncap + 4_i8b*nside*ir 856 else ! in South cap 857 npix = nside2npix(nside) 858 ir = 4_i8b*nside-ring - 1 ! count ring from south 859 npnorth = npix - ir*(ir+1_i8b)*2_i8b 860 endif 861 endif 862 return 863 end subroutine pixels_per_ring 864 !======================================================================= 865 subroutine check_edge_pixels(nside, nsboost, irmin, irmax, phi0, dphi, ringphi, ngr) 866 !======================================================================= 867 integer(i4b), intent(in) :: nside, nsboost, irmin, irmax 868 real(dp), intent(in) :: phi0 869 real(dp), dimension(1:), intent(in) :: dphi 870 integer(i4b), dimension(1:,1:), intent(inout) :: ringphi 871 integer(i4b), intent(inout) :: ngr 872 873 integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh 874 real(dp), dimension(1:2*nsboost+1) :: phiw, phie 875 real(dp) :: dd, dph, phic 876 877 !======================================================================= 878 if (nsboost <= 1) return 879 880 nrh = irmax-irmin+1 881 do i=1, ngr ! loop on low-res rings 882 i0 = ringphi(1,i) * nsboost - nsboost - irmin 883 do k=-1,1,2 ! West and East side of disc 884 kk = (k+5)/2 ! 2 or 3 885222 continue 886 iphi = ringphi(kk, i) 887 if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then 888 call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie) 889 do j=1, 2*nsboost+1 890 if (i0+j >= 1 .and. i0+j <= nrh) then 891 phic = (phie(j)+phiw(j))*0.5_dp ! pixel center 892 dph = (phie(j)-phiw(j))*0.5_dp + dphi(i0+j) ! pixel size + circle radius 893 dd = abs(phi0 - phic) ! distance from disc center to pixel border sample 894 dd = min(dd, twopi - dd) ! in [0,Pi] 895 if (dd <= dph) goto 1000 ! pixel touched by disc, move to next one 896 endif 897 enddo 898 ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards 899 goto 222 ! try next pixel inward 9001000 continue 901 endif 902 enddo ! loop on side 903 enddo ! loop on low-res rings 904 905 ! remove empty rings 906 ngr_out = 0 907 do i=1,ngr 908 diff = ringphi(3,i) - ringphi(2,i) 909 if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then 910 ngr_out = ngr_out + 1 911 ringphi(1:3, ngr_out) = ringphi(1:3, i) 912 endif 913 enddo 914 ! set empty rings to -1 915 do i=ngr_out+1, ngr 916 ringphi(2:3, i) = -1 917 enddo 918 919 ngr = ngr_out 920 return 921 end subroutine check_edge_pixels 922 !======================================================================= 923 subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie) 924 !======================================================================= 925 integer(i4b), intent(in) :: nside, nsboost, iring, iphi 926 real(dp), dimension(1:2*nsboost+1), intent(out) :: phiw, phie 927 928 real(dp), dimension(1:2*nsboost+1) :: f, f1, phiw_t, phie_t 929 real(dp) :: c0, quad, phie1, phie2, phiw1, phiw2, cv 930 integer(i4b) :: npr, kshift, nq, ip, i 931 logical(lgt) :: transition 932 !======================================================================= 933 934 call pixels_per_ring(nside, iring, npr, kshift) 935 !f = ((/ (i,i=0,2*nsboost) /) - nsboost) / nsboost 936 f = ((/ (i*1.d0,i=0,2*nsboost) /) - nsboost*1.d0) / nsboost 937 938 nq = npr/4 ! number of pixels on current ring in [0,Pi/2] (quadrant) 939 transition = (iring == nside .or. iring == nside*3) 940 941 if (nq == nside .or. transition) then ! equatorial region (and transition rings) 942 943 f1 = (1.0_dp-abs(f))*0.5_dp ! triangle of height 1/2 944 f1 = halfpi * f1 / nq 945 c0 = halfpi * (iphi + kshift*0.5_dp) / nq 946 phiw = c0 - f1 947 phie = c0 + f1 948 if (transition) then ! store for future use 949 phiw_t = phiw 950 phie_t = phie 951 endif 952 endif 953 954 if (nq < nside .or. transition) then ! polar regions and transition rings 955 ip = mod(iphi,nq) ! in [0,nq-1] 956 quad = iphi / nq ! quadrant in [0,3] 957 if (iring <= nside*2) then 958 f1 = halfpi / (nq + f) 959 else 960 f1 = halfpi / (nq - f)! swap sign for South pole 961 endif 962 do i=1, 2*nsboost+1 963 cv = f1(i) 964 phiw1 = min(cv * ip, halfpi) 965 phie1 = min(cv * ( ip+1), halfpi) 966 phiw2 = min(cv * (nq-ip-1), halfpi) 967 phie2 = min(cv * (nq-ip), halfpi) 968 phiw(i) = max(phiw1, halfpi - phie2) + (quad * halfpi) 969 phie(i) = min(phie1, halfpi - phiw2) + (quad * halfpi) 970 enddo 971 endif 972 973 if (transition) then 974 if (iring == nside) then ! transition in N hemisphere 975 phiw(nsboost+2:2*nsboost+1) = phiw_t(nsboost+2:2*nsboost+1) 976 phie(nsboost+2:2*nsboost+1) = phie_t(nsboost+2:2*nsboost+1) 977 else ! transition in S hemisphere 978 phiw(1:nsboost+1) = phiw_t(1:nsboost+1) 979 phie(1:nsboost+1) = phie_t(1:nsboost+1) 980 endif 981 endif 982 983 return 984 end subroutine find_pixel_bounds 985 986! !======================================================================= 987! subroutine correct_ring_phi(location, iring, iphi) 988! !======================================================================= 989! ! returns ring and phi indexes corrected from round-off errors 990! ! appearing at large Nside 991! ! if phi <0 : move 1 ring North 992! ! if phi >4*iring: move 1 ring South 993! ! rings are counted from the closest pole starting at 1. 994! integer(i4b), intent(in) :: location !+1:North, -1:South 995! integer(i4b), intent(inout) :: iring, iphi 996! integer(i4b) :: delta 997! !----------------------------------------------------------------------- 998! delta = 0 999! if (iphi < 0) delta = -1 1000! if (iphi >= 4*iring) delta = +1 1001! if (delta /= 0) then 1002! if (abs(location) /= 1) then 1003! stop 'wrong location' 1004! endif 1005! if (delta > 0) then ! too large iphi 1006! iphi = iphi - 4*iring ! use old iring 1007! iring = iring + location ! move south 1008! else ! too small iphi 1009! iring = iring - location ! move north 1010! iphi = iphi + 4*iring ! use new iring 1011! endif 1012! endif 1013! return 1014! end subroutine correct_ring_phi 1015 !======================================================================= 1016 ! CHEAP_ISQRT 1017 ! Returns exact Floor(sqrt(x)) where x is a (64 bit) integer. 1018 ! y^2 <= x < (y+1)^2 (1) 1019 ! The double precision floating point operation is not accurate enough 1020 ! when dealing with 64 bit integers, especially in the vicinity of 1021 ! perfect squares. 1022 !======================================================================= 1023#ifndef NO64BITS 1024 function cheap_isqrt_8(lin) result (lout) 1025 integer(i8b), intent(in) :: lin 1026 integer(i8b) :: lout, diff 1027 real(DP) :: dout, din 1028 lout = floor(sqrt(dble(lin)), kind=I8B) ! round-off error may offset result 1029 diff = lin - lout*lout ! test Eq (1) 1030 if (diff <0) lout = lout - 1 1031 if (diff >2*lout) lout = lout + 1 1032 return 1033 end function cheap_isqrt_8 1034#endif 1035 function cheap_isqrt_4(lin) result (lout) 1036 integer(i4b), intent(in) :: lin 1037 integer(i4b) :: lout 1038 lout = floor(sqrt(dble(lin)), kind=I4B) 1039 return 1040 end function cheap_isqrt_4 1041 !======================================================================= 1042 1043 ! perform 2 compilations of the same source file 1044#undef DOI8Bi 1045#include "pixel_routines.F90" 1046#ifndef NO64BITS 1047#define DOI8B 1048#include "pixel_routines.F90" 1049#endif 1050! 1051! pix2ang_ring 1052! pix2vec_ring 1053! ang2pix_ring 1054! vec2pix_ring 1055! pix2ang_nest 1056! pix2vec_nest 1057! ang2pix_nest 1058! vec2pix_nest 1059! nest2ring 1060! npix2nside 1061 1062 !************************************************************* 1063 ! 1064 ! MAP Manipulations 1065 ! 1066 !************************************************************* 1067 subroutine warning_oldbounds(code, cos_theta_cut, zbounds) 1068 character(len=*), intent(in) :: code 1069 real(DP), intent(in) :: cos_theta_cut 1070 real(DP), intent(out), dimension(1:2) :: zbounds 1071 1072 if (cos_theta_cut <= 0.0_dp) then ! no cut 1073 zbounds(1) = -1.0_dp 1074 zbounds(2) = 1.0_dp 1075 else 1076 zbounds(1) = cos_theta_cut 1077 zbounds(2) = - cos_theta_cut 1078 endif 1079 print*,' -------------------------------------' 1080 print*,'WARNING: obsolete interface to '//code 1081 print*,' cos_theta_cut currently a DP scalar with value' 1082 write(*,9000) ' ',cos_theta_cut 1083 print*,' shoud now be replaced with a 2-element vector with values:' 1084 write(*,9001) ' ',zbounds(1),zbounds(2) 1085 print*,' See documentation for details.' 1086 print*,' -------------------------------------' 10879000 format (a,g13.6) 10889001 format (a,g13.6,g13.6) 1089 1090 return 1091 end subroutine warning_oldbounds 1092 !============================================================== 1093 ! REMOVE_DIPOLE( nside, map, ordering, degree, multipoles, zbounds, fmissval, mask, weights) 1094 ! 1095 ! removes monopole (and dipole) from a map 1096 ! 1097 ! Nside: I4, IN : Healpix resolution parameter 1098 ! map: KMAP, array,INOUT: Heapix map (see Notes below) 1099 ! ordering: I4, IN: Healpix scheme 1:RING, 2: NESTED 1100 ! degree: I4, IN: multipole to remove, 1: monopole, 2: monopole and dipole 1101 ! multipoles:R8, array,OUT: value of monopole and dipole 1102 ! zbounds :R8, IN: 2-el vector 1103 ! range of z in [-1,1] on which to estimate monopole and dipole 1104 ! fmissval: KMAP, Option, IN: value used to flag bad pixel on input, default=-1.6375e30 1105 ! Pixels with map = fmissval are not used for fit 1106 ! mask : KMAP, Option, IN: Pixels with |mask|<1.e-10 are not used for fit 1107 ! others are kept as is 1108 ! note : the mask in NOT applied to the map 1109 ! weights : KMAP, Option, IN: pixel weight to be applied to the map 1110 ! 1111 ! KMAP: either R4 or R8 1112 ! 1113 ! note : if degree= 1, or 2, the map is modified on output 1114 ! * the monopole (and dipole) is/are removed 1115 ! * pixels within the symmetric cut parameterized 1116 ! by cos_theta_cut are set to fmissval (or its default value) 1117 ! if degree = 0, nothing is done 1118 ! all other values of degree are invalid 1119 ! 1120 ! v1.0, EH, Caltech, Jan-2002, based on homonyme IDL routine 1121 ! 2008-02-03: addition of weights 1122 ! 1123 !============================================================== 1124 subroutine remove_dipole_real( nside, map, ordering, degree, multipoles, zbounds, & 1125 & fmissval, mask, weights) 1126 !============================================================ 1127 use num_rec, only : dsvdcmp, dsvbksb 1128 ! parameters 1129 character(len=*), parameter :: code = "REMOVE_DIPOLE_REAL" 1130 integer(kind=I4B), parameter :: KMAP = SP 1131 real (kind=KMAP), parameter :: fbad_value = -1.6375e30_KMAP 1132 ! 1133 include 'remove_dipole_inc.f90' 1134 end subroutine remove_dipole_real 1135 !============================================================== 1136 subroutine remove_dipole_double( nside, map, ordering, degree, multipoles, zbounds, & 1137 & fmissval, mask, weights) 1138 !============================================================ 1139 use num_rec, only : dsvdcmp, dsvbksb 1140 ! parameters 1141 character(len=*), parameter :: code = "REMOVE_DIPOLE_DOUBLE" 1142 integer(kind=I4B), parameter :: KMAP = DP 1143 real (kind=KMAP), parameter :: fbad_value = -1.6375e30_KMAP 1144 ! 1145 include 'remove_dipole_inc.f90' 1146 end subroutine remove_dipole_double 1147 1148 !============================================================== 1149 subroutine remove_dipole_real_old( nside, map, ordering, degree, multipoles, cos_theta_cut, fmissval, mask) 1150 !============================================================ 1151 ! parameters 1152 character(len=*), parameter :: code = "REMOVE_DIPOLE_REAL" 1153 integer(kind=I4B), parameter :: KMAP = SP 1154 ! dummy 1155 integer(kind=i4b), intent(in) :: nside 1156 integer(kind=i4b), intent(in) :: ordering, degree 1157 real (kind=DP), dimension(0:), intent(out) :: multipoles 1158 real (kind=KMAP), dimension(0:), intent(inout) :: map 1159 real (kind=DP), intent(in) :: cos_theta_cut 1160 real (kind=KMAP), intent(in), optional :: fmissval 1161 real (kind=KMAP), dimension(0:), intent(in), optional :: mask 1162 ! local 1163 real(DP), dimension(1:2) :: zbounds 1164 1165 call warning_oldbounds(code, cos_theta_cut, zbounds) 1166 call remove_dipole(nside, map, ordering, degree, multipoles, zbounds, fmissval, mask) 1167 end subroutine remove_dipole_real_old 1168 !============================================================== 1169 subroutine remove_dipole_double_old( nside, map, ordering, degree, multipoles, cos_theta_cut, fmissval, mask) 1170 !============================================================ 1171 ! parameters 1172 character(len=*), parameter :: code = "REMOVE_DIPOLE_DOUBLE" 1173 integer(kind=I4B), parameter :: KMAP = DP 1174 ! dummy 1175 integer(kind=i4b), intent(in) :: nside 1176 integer(kind=i4b), intent(in) :: ordering, degree 1177 real (kind=DP), dimension(0:), intent(out) :: multipoles 1178 real (kind=KMAP), dimension(0:), intent(inout) :: map 1179 real (kind=DP), intent(in) :: cos_theta_cut 1180 real (kind=KMAP), intent(in), optional :: fmissval 1181 real (kind=KMAP), dimension(0:), intent(in), optional :: mask 1182 ! local 1183 real(DP), dimension(1:2) :: zbounds 1184 1185 call warning_oldbounds(code, cos_theta_cut, zbounds) 1186 call remove_dipole(nside, map, ordering, degree, multipoles, zbounds, fmissval, mask) 1187 end subroutine remove_dipole_double_old 1188 1189 !============================================================== 1190 subroutine remove_dipole_real_v12( nside, map, nmaps, ordering, degree, multipoles, cos_theta_cut, fmissval, mask) 1191 !============================================================ 1192 ! parameters 1193 character(len=*), parameter :: code = "REMOVE_DIPOLE_REAL" 1194 integer(kind=I4B), parameter :: KMAP = SP 1195 ! dummy 1196 integer(kind=i4b), intent(in) :: nside 1197 integer(kind=i4b), intent(in) :: ordering, degree, nmaps 1198 real (kind=DP), dimension(0:), intent(out) :: multipoles 1199 real (kind=KMAP), dimension(0:), intent(inout) :: map 1200 real (kind=DP), intent(in) :: cos_theta_cut 1201 real (kind=KMAP), intent(in), optional :: fmissval 1202 real (kind=KMAP), dimension(0:), intent(in), optional :: mask 1203 ! local 1204 real(DP), dimension(1:2) :: zbounds 1205 1206 print*,'==========================================================' 1207 print*,'WARNING: Interface to remove_dipole has changed' 1208 print*,' from' 1209 print*,'call remove_dipole(nside, map, NMAPS, ordering, degree, multipoles, COS_THETA_CUT, fmissval, mask)' 1210 print*,' to' 1211 print*,'call remove_dipole(nside, map, ordering, degree, multipoles, ZBOUNDS, fmissval, mask)' 1212 print*,'==========================================================' 1213 1214 call warning_oldbounds(code, cos_theta_cut, zbounds) 1215 call remove_dipole(nside, map, ordering, degree, multipoles, zbounds, fmissval, mask) 1216 end subroutine remove_dipole_real_v12 1217 !============================================================== 1218 subroutine remove_dipole_double_v12( nside, map, nmaps, ordering, degree, multipoles, cos_theta_cut, fmissval, mask) 1219 !============================================================ 1220 ! parameters 1221 character(len=*), parameter :: code = "REMOVE_DIPOLE_DOUBLE" 1222 integer(kind=I4B), parameter :: KMAP = DP 1223 ! dummy 1224 integer(kind=i4b), intent(in) :: nside 1225 integer(kind=i4b), intent(in) :: ordering, degree, nmaps 1226 real (kind=DP), dimension(0:), intent(out) :: multipoles 1227 real (kind=KMAP), dimension(0:), intent(inout) :: map 1228 real (kind=DP), intent(in) :: cos_theta_cut 1229 real (kind=KMAP), intent(in), optional :: fmissval 1230 real (kind=KMAP), dimension(0:), intent(in), optional :: mask 1231 ! local 1232 real(DP), dimension(1:2) :: zbounds 1233 1234 print*,'==========================================================' 1235 print*,'WARNING: Interface to remove_dipole has changed' 1236 print*,' from' 1237 print*,'call remove_dipole(nside, map, NMAPS, ordering, degree, multipoles, COS_THETA_CUT, fmissval, mask)' 1238 print*,' to' 1239 print*,'call remove_dipole(nside, map, ordering, degree, multipoles, ZBOUNDS, fmissval, mask)' 1240 print*,'==========================================================' 1241 1242 call warning_oldbounds(code, cos_theta_cut, zbounds) 1243 call remove_dipole(nside, map, ordering, degree, multipoles, zbounds, fmissval, mask) 1244 end subroutine remove_dipole_double_v12 1245 1246 !======================================================================= 1247 ! ADD_DIPOLE( nside, map, ordering, degree, multipoles, fmissval) 1248 ! 1249 ! removes monopole (and dipole) from a map 1250 ! 1251 ! Nside: I4, IN : Healpix resolution parameter 1252 ! map: KMAP, array,INOUT: Heapix map (see Notes below) 1253 ! ordering: I4, IN: Healpix scheme 1:RING, 2: NESTED 1254 ! degree: I4, IN: multipole to remove, 1: monopole, 2: monopole and dipole 1255 ! multipoles:R8, array,IN: value of monopole and dipole 1256 ! fmissval: KMAP, Option, IN: value used to flag bad pixel on input, default=-1.6375e30 1257 ! Pixels with map = fmissval are left unchanged 1258 ! 2009-03-25: accepts Nside > 8192 1259 !======================================================================= 1260 subroutine add_dipole_real(nside, map, ordering, degree, multipoles, fmissval) 1261 !======================================================================= 1262 ! single precision 1263 !======================================================================= 1264 character(len=*), parameter :: code = "ADD_DIPOLE_REAL" 1265 integer(kind=I4B), parameter :: KMAP = SP 1266 real (kind=KMAP), parameter :: fbad_value = -1.6375e30_KMAP 1267 ! dummy 1268 integer(kind=i4b), intent(in) :: nside 1269 integer(kind=i4b), intent(in) :: ordering, degree 1270 real (kind=DP), dimension(0:), intent(in) :: multipoles 1271 real (kind=KMAP), dimension(0:), intent(inout) :: map 1272 real (kind=KMAP), intent(in), optional :: fmissval 1273 1274 real (kind=KMAP) :: fmiss_effct 1275 integer(kind=i8b) :: ipix, npix 1276 logical(lgt) :: dodipole !, do_mask 1277 real(kind=dp), dimension(1:3) :: vec 1278 !======================================================================= 1279 1280 npix = nside2npix(nside) 1281 fmiss_effct = fbad_value 1282 if (present(fmissval)) fmiss_effct = fmissval 1283 1284 if (degree == 0) then 1285 print*," No monopole nor dipole to add" 1286 return 1287 elseif (degree == 1) then 1288 dodipole = .false. 1289 else if (degree == 2) then 1290 dodipole = .true. 1291 else 1292 print*,code//"> degree can only be " 1293 print*," 1: monopole (l=0) addition or " 1294 print*," 2: monopole and dipole (l=0,1) addition" 1295 print*,code//"> ABORT ! " 1296 call fatal_error 1297 endif 1298 1299 do ipix = 0, npix-1 1300 if ( abs(map(ipix) - fmiss_effct) <= abs(1.e-5*fmiss_effct) ) goto 20 1301! if (do_mask) then 1302! if (abs(mask(ipix)) <= 1.e-10) goto 20 1303! endif 1304 map(ipix) = map(ipix) + multipoles(0) 1305 if (dodipole) then 1306 ! computes dipole basis functions 1307 ! pixel -> vector 1308 if (ordering == 1) call pix2vec_ring( nside, ipix, vec) 1309 if (ordering == 2) call pix2vec_nest( nside, ipix, vec) 1310 map(ipix) = map(ipix) + sum(multipoles(1:3) * vec(1:3)) 1311 endif 1312 131320 continue 1314 enddo 1315 1316 return 1317 end subroutine add_dipole_real 1318 !======================================================================= 1319 subroutine add_dipole_double(nside, map, ordering, degree, multipoles, fmissval) 1320 !======================================================================= 1321 ! single precision 1322 !======================================================================= 1323 character(len=*), parameter :: code = "ADD_DIPOLE_DOUBLE" 1324 integer(kind=I4B), parameter :: KMAP = DP 1325 real (kind=KMAP), parameter :: fbad_value = -1.6375e30_KMAP 1326 ! dummy 1327 integer(kind=i4b), intent(in) :: nside 1328 integer(kind=i4b), intent(in) :: ordering, degree 1329 real (kind=DP), dimension(0:), intent(in) :: multipoles 1330 real (kind=KMAP), dimension(0:), intent(inout) :: map 1331 real (kind=KMAP), intent(in), optional :: fmissval 1332 1333 real (kind=KMAP) :: fmiss_effct 1334 integer(kind=i8b) :: ipix, npix 1335 logical(lgt) :: dodipole !, do_mask 1336 real(kind=dp), dimension(1:3) :: vec 1337 !======================================================================= 1338 1339 npix = nside2npix(nside) 1340 fmiss_effct = fbad_value 1341 if (present(fmissval)) fmiss_effct = fmissval 1342 1343 if (degree == 0) then 1344 print*," No monopole nor dipole to add" 1345 return 1346 elseif (degree == 1) then 1347 dodipole = .false. 1348 else if (degree == 2) then 1349 dodipole = .true. 1350 else 1351 print*,code//"> degree can only be " 1352 print*," 1: monopole (l=0) addition or " 1353 print*," 2: monopole and dipole (l=0,1) addition" 1354 print*,code//"> ABORT ! " 1355 call fatal_error 1356 endif 1357 1358 do ipix = 0, npix-1 1359 if ( abs(map(ipix) - fmiss_effct) <= abs(1.e-5*fmiss_effct) ) goto 20 1360! if (do_mask) then 1361! if (abs(mask(ipix)) <= 1.e-10) goto 20 1362! endif 1363 map(ipix) = map(ipix) + multipoles(0) 1364 if (dodipole) then 1365 ! computes dipole basis functions 1366 ! pixel -> vector 1367 if (ordering == 1) call pix2vec_ring( nside, ipix, vec) 1368 if (ordering == 2) call pix2vec_nest( nside, ipix, vec) 1369 map(ipix) = map(ipix) + sum(multipoles(1:3) * vec(1:3)) 1370 endif 1371 137220 continue 1373 enddo 1374 1375 return 1376 end subroutine add_dipole_double 1377 1378 !==================================================== 1379 ! apply_mask (map, order, mask=mask, zbounds=zbounds) 1380 ! applies a mask and/or a zbounds to a TQU maps, 1381 ! map and mask should have the same ordering 1382 ! map : SP/DP, INOUT, (0:,1:) 1383 ! order: I4B, IN, 1=RING, 2=NESTED 1384 ! mask: SP/DP IN, (0:,1:) optional 1385 ! zbounds: DP IN, (1:2) optional 1386 !===================================================================== 1387 subroutine apply_mask_real(map, order, mask, zbounds) 1388 !===================================================================== 1389 integer(I4B), parameter :: KMAP = SP 1390 include 'apply_mask_inc.f90' 1391 end subroutine apply_mask_real 1392 !===================================================================== 1393 subroutine apply_mask_double(map, order, mask, zbounds) 1394 !===================================================================== 1395 integer(I4B), parameter :: KMAP = DP 1396 include 'apply_mask_inc.f90' 1397 end subroutine apply_mask_double 1398 !==================================================== 1399 ! medfiltmap 1400 ! compute the median filtered map of a given Healpix map 1401 ! in_map: SP/DP input Healpix full sky map 1402 ! radius: DP radius in Radians 1403 ! med_map: SP/DP output Healpix full sky map 1404 ! nest: I4B, OPT either 0 (ring scheme) or 1 (nested scheme) 1405 ! fmissval:SP/DP, OPT sentinel value given to missing pixels 1406 ! fill_holes: LGT, OPT 1407 ! 2012-07-17: Parallel OpenMP implementation 1408 !================================================================= 1409 subroutine medfiltmap_S( in_map, radius, med_map, nest, fmissval, fill_holes) 1410 !================================================================= 1411 use statistics, only: median 1412 integer(I4B), parameter :: KMAP = SP 1413 ! 1414 real(KMAP), dimension(0:),intent(in), target :: in_map 1415 real(DP), intent(in) :: radius 1416 real(KMAP), dimension(0:),intent(out) :: med_map 1417 integer(I4B), intent(in), optional :: nest 1418 real(KMAP), intent(in), optional :: fmissval 1419 logical(LGT), intent(in), optional :: fill_holes 1420 ! 1421 integer(I8B) :: npix, np 1422 integer(I4B) :: nside, p, nlist, status 1423 logical(LGT) :: do_nest, do_fill 1424 real(DP) :: fraction 1425 real(KMAP) :: fmissval_in 1426 integer(I4B), dimension(:), allocatable :: listpix 1427 real(DP), dimension(1:3) :: vector 1428 character(len=*), parameter :: code = 'medfiltmap' 1429 !----------------------------------------------- 1430 npix = long_size(in_map) 1431 nside = npix2nside(npix) 1432 call assert(nside > 0, code//": invalid map size") 1433 1434 fraction = 0.5_DP * (1.0_dp - cos(radius)) 1435 np = npix * fraction * 1.2 + 50 1436 call assert(np < MAX_I4B, code//": too many pixels to compute median") 1437 1438 do_nest = .false. 1439 if (present(nest)) then 1440 call assert(nest>=0 .and. nest <=1,code//': invalid NEST flag') 1441 do_nest = (nest == 1) 1442 endif 1443 1444 do_fill = .false. 1445 if (present(fill_holes)) do_fill = fill_holes 1446 1447 fmissval_in = hpx_Sbadval 1448 if (present(fmissval)) fmissval_in = fmissval 1449 1450 ! make sure common arrays are initiated 1451 call mk_pix2xy() 1452!!! print*,'************* Parallel Median **************' 1453!$OMP parallel default(none) & 1454!$OMP shared(in_map, med_map, pix2x, pix2y, & 1455!$OMP nside, npix, radius, np, do_nest, do_fill, nest, fmissval_in) & 1456!$OMP private(listpix, vector, p, nlist, status) 1457 allocate(listpix(0:np-1),stat=status) 1458 call assert_alloc(status,code,'listpix') 1459!$OMP do schedule(dynamic, 64) 1460 do p = 0, npix-1 1461 ! find pixel location 1462 if (do_nest) then 1463 call pix2vec_nest( nside, p, vector) 1464 else 1465 call pix2vec_ring( nside, p, vector) 1466 endif 1467 1468 ! find disc centered on pixel 1469 call query_disc(nside, vector, radius, listpix, nlist, nest=nest) 1470 1471 if (do_fill .or. abs(in_map(p)-fmissval_in) > abs(fmissval_in*1.e-7)) then 1472 med_map(p) = median(in_map(listpix(0:nlist-1)), badval = fmissval_in, even=(nlist<100)) 1473 else 1474 med_map(p) = in_map(p) 1475 endif 1476 enddo 1477!$OMP end do 1478 deallocate(listpix) 1479!$OMP end parallel 1480 return 1481 end subroutine medfiltmap_S 1482 !================================================================= 1483 subroutine medfiltmap_D( in_map, radius, med_map, nest, fmissval, fill_holes) 1484 !================================================================= 1485 use statistics, only: median 1486 integer(I4B), parameter :: KMAP = DP 1487 ! 1488 real(KMAP), dimension(0:),intent(in), target :: in_map 1489 real(DP), intent(in) :: radius 1490 real(KMAP), dimension(0:),intent(out) :: med_map 1491 integer(I4B), intent(in), optional :: nest 1492 real(KMAP), intent(in), optional :: fmissval 1493 logical(LGT), intent(in), optional :: fill_holes 1494 ! 1495 integer(I8B) :: npix, np 1496 integer(I4B) :: nside, p, nlist, status 1497 logical(LGT) :: do_nest, do_fill 1498 real(DP) :: fraction 1499 real(KMAP) :: fmissval_in 1500 integer(I4B), dimension(:), allocatable :: listpix 1501 real(DP), dimension(1:3) :: vector 1502 character(len=*), parameter :: code = 'medfiltmap' 1503 !----------------------------------------------- 1504 npix = long_size(in_map) 1505 nside = npix2nside(npix) 1506 call assert(nside > 0, code//": invalid map size") 1507 1508 fraction = 0.5_DP * (1.0_dp - cos(radius)) 1509 np = npix * fraction * 1.2 + 50 1510 call assert(np < MAX_I4B, code//": too many pixels to compute median") 1511 1512 do_nest = .false. 1513 if (present(nest)) then 1514 call assert(nest>=0 .and. nest <=1,code//': invalid NEST flag') 1515 do_nest = (nest == 1) 1516 endif 1517 1518 do_fill = .false. 1519 if (present(fill_holes)) do_fill = fill_holes 1520 1521 fmissval_in = hpx_Dbadval 1522 if (present(fmissval)) fmissval_in = fmissval 1523 1524 ! make sure common arrays are initiated 1525 call mk_pix2xy() 1526!!! print*,'************* Parallel Median **************' 1527!$OMP parallel default(none) & 1528!$OMP shared(in_map, med_map, pix2x, pix2y, & 1529!$OMP nside, npix, radius, np, do_nest, do_fill, nest, fmissval_in) & 1530!$OMP private(listpix, vector, p, nlist, status) 1531 allocate(listpix(0:np-1),stat=status) 1532 call assert_alloc(status,code,'listpix') 1533!$OMP do schedule(dynamic, 64) 1534 do p = 0, npix-1 1535 ! find pixel location 1536 if (do_nest) then 1537 call pix2vec_nest( nside, p, vector) 1538 else 1539 call pix2vec_ring( nside, p, vector) 1540 endif 1541 1542 ! find disc centered on pixel 1543 call query_disc(nside, vector, radius, listpix, nlist, nest=nest) 1544 1545 if (do_fill .or. abs(in_map(p)-fmissval_in) > abs(fmissval_in*1.e-7)) then 1546 med_map(p) = median(in_map(listpix(0:nlist-1)), badval = fmissval_in, even=(nlist<100)) 1547 else 1548 med_map(p) = in_map(p) 1549 endif 1550 enddo 1551!$OMP end do 1552 deallocate(listpix) 1553!$OMP end parallel 1554 return 1555 end subroutine medfiltmap_D 1556 !************************************************************** 1557 ! following 2x3 routines: out of place conversions for 1D maps 1558 ! if size(map) = Npix, peak memory = 2*Npix 1559 ! These routines are parallelized for shared memory architecture 1560 ! (using OpenMP directives) 1561 !************************************************************** 1562 !======================================================================= 1563 ! makes the conversion NEST to RING 1564 !======================================================================= 1565 subroutine convert_nest2ring_int_1d(nside, map) 1566 !======================================================================= 1567 character(len=*), parameter :: code = 'convert_nest2ring_int_1d' 1568 integer(kind=I4B), parameter :: KMAP = I4B 1569 integer(kind=KMAP), dimension(0:), intent(inout) :: map 1570 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1571 include 'convert_nest2ring_1d_inc.f90' 1572 end subroutine convert_nest2ring_int_1d 1573#ifndef NO64BITS 1574 !======================================================================= 1575 subroutine convert_nest2ring_int8_1d(nside, map) 1576 !======================================================================= 1577 character(len=*), parameter :: code = 'convert_nest2ring_int8_1d' 1578 integer(kind=I4B), parameter :: KMAP = I8B 1579 integer(kind=KMAP), dimension(0:), intent(inout) :: map 1580 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1581 include 'convert_nest2ring_1d_inc.f90' 1582 end subroutine convert_nest2ring_int8_1d 1583#endif 1584 !======================================================================= 1585 subroutine convert_nest2ring_real_1d(nside, map) 1586 !======================================================================= 1587 character(len=*), parameter :: code = 'convert_nest2ring_real_1d' 1588 integer(kind=I4B), parameter :: KMAP = SP 1589 real (kind=KMAP), dimension(0:), intent(inout) :: map 1590 real (kind=KMAP), dimension(:), allocatable :: map_tmp 1591 include 'convert_nest2ring_1d_inc.f90' 1592 end subroutine convert_nest2ring_real_1d 1593 !======================================================================= 1594 subroutine convert_nest2ring_double_1d(nside, map) 1595 !======================================================================= 1596 character(len=*), parameter :: code = 'convert_nest2ring_double_1d' 1597 integer(kind=I4B), parameter :: KMAP = DP 1598 real (kind=KMAP), dimension(0:), intent(inout) :: map 1599 real (kind=KMAP), dimension(:), allocatable :: map_tmp 1600 include 'convert_nest2ring_1d_inc.f90' 1601 end subroutine convert_nest2ring_double_1d 1602 !======================================================================= 1603 ! makes the conversion RING to NEST 1604 !======================================================================= 1605 subroutine convert_ring2nest_int_1d(nside, map) 1606 !======================================================================= 1607 character(len=*), parameter :: code = 'convert_ring2nest_int_1d' 1608 integer(kind=I4B), parameter :: KMAP = I4B 1609 integer(kind=KMAP), dimension(0:), intent(inout) :: map 1610 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1611 include 'convert_ring2nest_1d_inc.f90' 1612 end subroutine convert_ring2nest_int_1d 1613#ifndef NO64BITS 1614 !======================================================================= 1615 subroutine convert_ring2nest_int8_1d(nside, map) 1616 !======================================================================= 1617 character(len=*), parameter :: code = 'convert_ring2nest_int8_1d' 1618 integer(kind=I4B), parameter :: KMAP = I8B 1619 integer(kind=KMAP), dimension(0:), intent(inout) :: map 1620 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1621 include 'convert_ring2nest_1d_inc.f90' 1622 end subroutine convert_ring2nest_int8_1d 1623#endif 1624 !======================================================================= 1625 subroutine convert_ring2nest_real_1d(nside, map) 1626 !======================================================================= 1627 character(len=*), parameter :: code = 'convert_ring2nest_real_1d' 1628 integer(kind=I4B), parameter :: KMAP = SP 1629 real (kind=KMAP), dimension(0:), intent(inout) :: map 1630 real (kind=KMAP), dimension(:), allocatable :: map_tmp 1631 include 'convert_ring2nest_1d_inc.f90' 1632 end subroutine convert_ring2nest_real_1d 1633 !======================================================================= 1634 subroutine convert_ring2nest_double_1d(nside, map) 1635 !======================================================================= 1636 character(len=*), parameter :: code = 'convert_ring2nest_double_1d' 1637 integer(kind=I4B), parameter :: KMAP = DP 1638 real (kind=KMAP), dimension(0:), intent(inout) :: map 1639 real (kind=KMAP), dimension(:), allocatable :: map_tmp 1640 include 'convert_ring2nest_1d_inc.f90' 1641 end subroutine convert_ring2nest_double_1d 1642 1643 !************************************************************** 1644 ! following 6 routines: out of place conversions for N-Dim maps 1645 ! if size(map) = Npix*Nd, peak memory = (Nd+2)*Npix 1646 ! 2004-08-25, EH 1647 !************************************************************** 1648 !======================================================================= 1649 subroutine convert_nest2ring_int_nd(nside, map) 1650 !======================================================================= 1651 ! NEST to RING conversion: 1D, integer 1652 !======================================================================= 1653 character(len=*), parameter :: code = "convert_nest2ring_int_nd" 1654 integer(kind=I4B), parameter :: KMAP = I4B 1655 integer(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1656 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1657 integer(kind=KMAP), dimension(:), pointer :: map1 1658 include 'convert_nest2ring_nd_inc.f90' 1659 end subroutine convert_nest2ring_int_nd 1660#ifndef NO64BITS 1661 !======================================================================= 1662 subroutine convert_nest2ring_int8_nd(nside, map) 1663 !======================================================================= 1664 ! NEST to RING conversion: 1D, integer*8 1665 !======================================================================= 1666 character(len=*), parameter :: code = "convert_nest2ring_int_nd" 1667 integer(kind=I4B), parameter :: KMAP = I8B 1668 integer(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1669 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1670 integer(kind=KMAP), dimension(:), pointer :: map1 1671 include 'convert_nest2ring_nd_inc.f90' 1672 end subroutine convert_nest2ring_int8_nd 1673#endif 1674 !======================================================================= 1675 subroutine convert_nest2ring_real_nd(nside, map) 1676 !======================================================================= 1677 ! NEST to RING conversion: 1D, real 1678 !======================================================================= 1679 character(len=*), parameter :: code = "convert_nest2ring_real_nd" 1680 integer(kind=I4B), parameter :: KMAP = SP 1681 real(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1682 real(kind=KMAP), dimension(:), allocatable :: map_tmp 1683 real(kind=KMAP), dimension(:), pointer :: map1 1684 include 'convert_nest2ring_nd_inc.f90' 1685 end subroutine convert_nest2ring_real_nd 1686 !======================================================================= 1687 subroutine convert_nest2ring_double_nd(nside, map) 1688 !======================================================================= 1689 ! NEST to RING conversion: 1D, double 1690 !======================================================================= 1691 character(len=*), parameter :: code = "convert_nest2ring_double_nd" 1692 integer(kind=I4B), parameter :: KMAP = DP 1693 real(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1694 real(kind=KMAP), dimension(:), allocatable :: map_tmp 1695 real(kind=KMAP), dimension(:), pointer :: map1 1696 include 'convert_nest2ring_nd_inc.f90' 1697 end subroutine convert_nest2ring_double_nd 1698 !======================================================================= 1699 subroutine convert_ring2nest_int_nd(nside, map) 1700 !======================================================================= 1701 ! RING to NEST conversion: 1D, integer 1702 !======================================================================= 1703 character(len=*), parameter :: code = "convert_ring2nest_int_nd" 1704 integer(kind=I4B), parameter :: KMAP = I4B 1705 integer(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1706 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1707 integer(kind=KMAP), dimension(:), pointer :: map1 1708 include 'convert_ring2nest_nd_inc.f90' 1709 end subroutine convert_ring2nest_int_nd 1710#ifndef NO64BITS 1711 !======================================================================= 1712 subroutine convert_ring2nest_int8_nd(nside, map) 1713 !======================================================================= 1714 ! RING to NEST conversion: 1D, integer 1715 !======================================================================= 1716 character(len=*), parameter :: code = "convert_ring2nest_int8_nd" 1717 integer(kind=I4B), parameter :: KMAP = I8B 1718 integer(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1719 integer(kind=KMAP), dimension(:), allocatable :: map_tmp 1720 integer(kind=KMAP), dimension(:), pointer :: map1 1721 include 'convert_ring2nest_nd_inc.f90' 1722 end subroutine convert_ring2nest_int8_nd 1723#endif 1724 !======================================================================= 1725 subroutine convert_ring2nest_real_nd(nside, map) 1726 !======================================================================= 1727 ! RING to NEST conversion: 1D, real 1728 !======================================================================= 1729 character(len=*), parameter :: code = "convert_ring2nest_real_nd" 1730 integer(kind=I4B), parameter :: KMAP = SP 1731 real(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1732 real(kind=KMAP), dimension(:), allocatable :: map_tmp 1733 real(kind=KMAP), dimension(:), pointer :: map1 1734 include 'convert_ring2nest_nd_inc.f90' 1735 end subroutine convert_ring2nest_real_nd 1736 !======================================================================= 1737 subroutine convert_ring2nest_double_nd(nside, map) 1738 !======================================================================= 1739 ! RING to NEST conversion: 1D, double 1740 !======================================================================= 1741 character(len=*), parameter :: code = "convert_ring2nest_double_nd" 1742 integer(kind=I4B), parameter :: KMAP = DP 1743 real(kind=KMAP), dimension(0:,1:), intent(inout), target :: map 1744 real(kind=KMAP), dimension(:), allocatable :: map_tmp 1745 real(kind=KMAP), dimension(:), pointer :: map1 1746 include 'convert_ring2nest_nd_inc.f90' 1747 end subroutine convert_ring2nest_double_nd 1748 !==================================================================== 1749 ! The following 6 routines convert in place integer, real, and double 1750 ! arrays between the NESTED and RING schemes. 1751 ! 1752 ! in place: without allocating a temporary map. This routine is more general, 1753 ! but slower than convert_nest2ring. 1754 ! 1755 ! This is a wrapper for the toolbox functions "ring2nest" and 1756 ! "nest2ring". Their names are supplied in the "subcall" 1757 ! argument. 1758 ! 1759 ! Author: Benjamin D. Wandelt October 1997 1760 ! Added to pix_tools for version 1.00 in March 1999. 1761 ! 2004-08-25: EH, added N-Dim facility and double precision IO 1762 !==================================================================== 1763 !==================================================================== 1764 subroutine convert_inplace_int_1d(subcall,map) 1765 !================================================================== 1766 ! 1D, integer implementation 1767 !================================================================== 1768 character(len=*), parameter :: code = "convert_inplace_int_1d" 1769 integer(kind=I4B), parameter :: KMAP = I4B 1770 integer(kind=KMAP), dimension(0:) :: map 1771 integer(kind=KMAP) :: pixbuf1,pixbuf2 1772 include 'convert_inplace_1d_inc.f90' 1773 end subroutine convert_inplace_int_1d 1774 !==================================================================== 1775 subroutine convert_inplace_real_1d(subcall,map) 1776 !==================================================================== 1777 ! 1D, real implementation 1778 !================================================================== 1779 character(len=*), parameter :: code = "convert_inplace_real_1d" 1780 integer(kind=I4B), parameter :: KMAP = SP 1781 real (kind=KMAP), dimension(0:) :: map 1782 real (kind=KMAP) :: pixbuf1,pixbuf2 1783 include 'convert_inplace_1d_inc.f90' 1784 end subroutine convert_inplace_real_1d 1785 !==================================================================== 1786 subroutine convert_inplace_double_1d(subcall,map) 1787 !==================================================================== 1788 ! 1D, double precision implementation 1789 !================================================================== 1790 character(len=*), parameter :: code = "convert_inplace_double_1d" 1791 integer(kind=I4B), parameter :: KMAP = DP 1792 real (kind=KMAP), dimension(0:) :: map 1793 real (kind=KMAP) :: pixbuf1,pixbuf2 1794 include 'convert_inplace_1d_inc.f90' 1795 end subroutine convert_inplace_double_1d 1796 !==================================================================== 1797 subroutine convert_inplace_int_nd(subcall,map) 1798 !================================================================== 1799 ! ND, integer implementation 1800 !================================================================== 1801 character(len=*), parameter :: code = "convert_inplace_int_nd" 1802 integer(kind=i4b), parameter :: ND_MAX = 10 1803 integer(kind=I4B), parameter :: KMAP = I4B 1804 integer(kind=KMAP), dimension(0:,1:) :: map 1805 integer(kind=KMAP), dimension(1:ND_MAX) :: pixbuf1,pixbuf2 1806 include 'convert_inplace_nd_inc.f90' 1807 end subroutine convert_inplace_int_nd 1808 !==================================================================== 1809 subroutine convert_inplace_real_nd(subcall,map) 1810 !==================================================================== 1811 ! ND, real implementation 1812 !================================================================== 1813 character(len=*), parameter :: code = "convert_inplace_real_nd" 1814 integer(kind=i4b), parameter :: ND_MAX = 10 1815 integer(kind=I4B), parameter :: KMAP = SP 1816 real (kind=KMAP), dimension(0:,1:) :: map 1817 real (kind=KMAP), dimension(1:ND_MAX) :: pixbuf1,pixbuf2 1818 include 'convert_inplace_nd_inc.f90' 1819 end subroutine convert_inplace_real_nd 1820 !==================================================================== 1821 subroutine convert_inplace_double_nd(subcall,map) 1822 !==================================================================== 1823 ! ND, double precision implementation 1824 !================================================================== 1825 character(len=*), parameter :: code = "convert_inplace_double_nd" 1826 integer(kind=i4b), parameter :: ND_MAX = 10 1827 integer(kind=I4B), parameter :: KMAP = DP 1828 real (kind=KMAP), dimension(0:,1:) :: map 1829 real (kind=KMAP), dimension(1:ND_MAX) :: pixbuf1,pixbuf2 1830 include 'convert_inplace_nd_inc.f90' 1831 end subroutine convert_inplace_double_nd 1832 !======================================================================= 1833 subroutine mk_pix2xy() 1834 !======================================================================= 1835 ! constructs the array giving x and y in the face from pixel number 1836 ! for the nested (quad-cube like) ordering of pixels 1837 ! 1838 ! the bits corresponding to x and y are interleaved in the pixel number 1839 ! one breaks up the pixel number by even and odd bits 1840 !======================================================================= 1841 INTEGER(KIND=I4B) :: kpix, jpix, ix, iy, ip, id 1842 1843 !cc cf block data data pix2x(1023) /0/ 1844 !----------------------------------------------------------------------- 1845 ! print *, 'initiate pix2xy' 1846 do kpix=0,1023 ! pixel number 1847 jpix = kpix 1848 IX = 0 1849 IY = 0 1850 IP = 1 ! bit position (in x and y) 1851! do while (jpix/=0) ! go through all the bits 1852 do 1853 if (jpix == 0) exit ! go through all the bits 1854 ID = MODULO(jpix,2) ! bit value (in kpix), goes in ix 1855 jpix = jpix/2 1856 IX = ID*IP+IX 1857 1858 ID = MODULO(jpix,2) ! bit value (in kpix), goes in iy 1859 jpix = jpix/2 1860 IY = ID*IP+IY 1861 1862 IP = 2*IP ! next bit (in x and y) 1863 enddo 1864 pix2x(kpix) = IX ! in 0,31 1865 pix2y(kpix) = IY ! in 0,31 1866 enddo 1867 1868 return 1869 end subroutine mk_pix2xy 1870! !======================================================================= 1871! subroutine mk_xy2pix() 1872! !======================================================================= 1873! ! sets the array giving the number of the pixel lying in (x,y) 1874! ! x and y are in {1,128} 1875! ! the pixel number is in {0,128**2-1} 1876! ! 1877! ! if i-1 = sum_p=0 b_p * 2^p 1878! ! then ix = sum_p=0 b_p * 4^p 1879! ! iy = 2*ix 1880! ! ix + iy in {0, 128**2 -1} 1881! !======================================================================= 1882! INTEGER(KIND=I4B):: k,ip,i,j,id 1883! !======================================================================= 1884 1885! do i = 1,128 !for converting x,y into 1886! j = i-1 !pixel numbers 1887! k = 0 1888! ip = 1 1889 1890! do 1891! if (j==0) then 1892! x2pix(i) = k 1893! y2pix(i) = 2*k 1894! exit 1895! else 1896! id = MODULO(J,2) 1897! j = j/2 1898! k = ip*id+k 1899! ip = ip*4 1900! endif 1901! enddo 1902 1903! enddo 1904 1905! return 1906! end subroutine mk_xy2pix 1907 !======================================================================= 1908 subroutine mk_xy2pix1() 1909 !======================================================================= 1910 ! sets the array giving the number of the pixel lying in (x,y) 1911 ! x and y are in {1,128} 1912 ! the pixel number is in {0,128**2-1} 1913 ! 1914 ! if i-1 = sum_p=0 b_p * 2^p 1915 ! then ix = sum_p=0 b_p * 4^p 1916 ! iy = 2*ix 1917 ! ix + iy in {0, 128**2 -1} 1918 !======================================================================= 1919 INTEGER(KIND=I4B):: k,ip,i,j,id 1920 !======================================================================= 1921 1922 do i = 0,127 !for converting x,y into 1923 j = i !pixel numbers 1924 k = 0 1925 ip = 1 1926 1927 do 1928 if (j==0) then 1929 x2pix1(i) = k 1930 y2pix1(i) = 2*k 1931 exit 1932 else 1933 id = MODULO(J,2) 1934 j = j/2 1935 k = ip*id+k 1936 ip = ip*4 1937 endif 1938 enddo 1939 1940 enddo 1941 1942 return 1943 end subroutine mk_xy2pix1 1944 !======================================================================= 1945 1946 !======================================================================= 1947 subroutine ang2vec(theta, phi, vector) 1948 !======================================================================= 1949 ! renders the vector (x,y,z) corresponding to angles 1950 ! theta (co-latitude measured from North pole, in [0,Pi] radians) 1951 ! and phi (longitude measured eastward, in radians) 1952 ! North pole is (x,y,z)=(0,0,1) 1953 ! added by EH, Feb 2000 1954 !======================================================================= 1955 REAL(KIND=DP), INTENT(IN) :: theta, phi 1956 REAL(KIND=DP), INTENT(OUT), dimension(1:) :: vector 1957 1958 REAL(KIND=DP) :: sintheta 1959 !======================================================================= 1960 1961 if (theta<0.0_dp .or. theta>pi) then 1962 print*,"ANG2VEC: theta : ",theta," is out of range [0, Pi]" 1963 call fatal_error 1964 endif 1965 sintheta = SIN(theta) 1966 1967 vector(1) = sintheta * COS(phi) 1968 vector(2) = sintheta * SIN(phi) 1969 vector(3) = COS(theta) 1970 1971 return 1972 end subroutine ang2vec 1973 !======================================================================= 1974 subroutine vec2ang(vector, theta, phi) 1975 !======================================================================= 1976 ! renders the angles theta, phi corresponding to vector (x,y,z) 1977 ! theta (co-latitude measured from North pole, in [0,Pi] radians) 1978 ! and phi (longitude measured eastward, in [0,2Pi[ radians) 1979 ! North pole is (x,y,z)=(0,0,1) 1980 ! added by EH, Feb 2000 1981 ! 2011-08: replaced ACOS(z) by more accurate ATAN2(r,z) 1982 !======================================================================= 1983 REAL(KIND=DP), INTENT(IN), dimension(1:) :: vector 1984 REAL(KIND=DP), INTENT(OUT) :: theta, phi 1985 1986 !======================================================================= 1987 1988 theta = atan2(sqrt(vector(1)**2 + vector(2)**2), vector(3)) 1989 1990 phi = 0.0_dp 1991 if (vector(1) /= 0.0_dp .or. vector(2) /= 0.0_dp) & 1992 & phi = ATAN2(vector(2),vector(1)) ! phi in ]-pi,pi] 1993 if (phi < 0.0) phi = phi + twopi ! phi in [0,2pi[ 1994 1995 return 1996 end subroutine vec2ang 1997 !======================================================================= 1998 function nside2npix(nside) result(npix_result) 1999 !======================================================================= 2000 ! given nside, returns npix such that npix = 12*nside^2 2001 ! nside should be a power of 2 smaller than ns_max 2002 ! if not, -1 is returned 2003 ! EH, Feb-2000 2004 ! 2009-03-04: returns i8b result, faster 2005 !======================================================================= 2006 INTEGER(KIND=I8B) :: npix_result 2007 INTEGER(KIND=I4B), INTENT(IN) :: nside 2008 2009 INTEGER(KIND=I8B) :: npix 2010 CHARACTER(LEN=*), PARAMETER :: code = "nside2npix" 2011 !======================================================================= 2012 2013 npix = (12_i8b*nside)*nside 2014 if (nside < 1 .or. nside > ns_max .or. iand(nside-1,nside) /= 0) then 2015 print*,code//": Nside="//trim(string(nside))//" is not a power of 2." 2016 npix = -1 2017 endif 2018 npix_result = npix 2019 2020 return 2021 end function nside2npix 2022 !======================================================================= 2023 subroutine surface_triangle(vec1, vec2, vec3, surface) 2024 !======================================================================= 2025 ! returns the surface in steradians 2026 ! of the spherical triangle with vertices vec1, vec2, vec3 2027 ! 2028 ! algorithm : finds triangle sides and uses l'Huilier formula to compute 2029 ! "spherical excess" = surface area of triangle on a sphere of radius one 2030 ! see, eg Bronshtein, Semendyayev Eq 2.86 2031 !======================================================================= 2032 real(kind=dp), dimension(1:), intent(in) :: vec1, vec2, vec3 2033 real(kind=dp), intent(out) :: surface 2034 2035 ! real(kind=dp), dimension(1:3) :: v1, v2, v3 2036 real(kind=dp), dimension(1:3) :: side 2037 ! real(kind=dp) :: hp 2038 real(kind=dp) :: x0, x1, x2, x3 2039 !======================================================================= 2040 2041 ! half perimeter 2042! hp = 0.5_dp * (side1 + side2 + side3) 2043 2044! ! l'Huilier formula 2045! x0 = tan( hp * 0.5_dp) 2046! x1 = tan((hp - side1) * 0.5_dp) 2047! x2 = tan((hp - side2) * 0.5_dp) 2048! x3 = tan((hp - side3) * 0.5_dp) 2049 2050 ! find triangle sides 2051 call angdist(vec2, vec3, side(1)) 2052 call angdist(vec3, vec1, side(2)) 2053 call angdist(vec1, vec2, side(3)) 2054 ! divide by 4 2055 side(1:3) = side(1:3) * 0.25_dp 2056 2057 ! l'Huilier formula 2058 x0 = tan( side(1) + side(2) + side(3) ) 2059 x1 = tan(-side(1) + side(2) + side(3) ) 2060 x2 = tan( side(1) - side(2) + side(3) ) 2061 x3 = tan( side(1) + side(2) - side(3) ) 2062 surface = 4.0_dp * atan( sqrt(x0 * x1 * x2 * x3) ) 2063 2064 return 2065 end subroutine surface_triangle 2066 !======================================================================= 2067 subroutine angdist(v1, v2, dist) 2068 !======================================================================= 2069 ! call angdist(v1, v2, dist) 2070 ! computes the angular distance dist (in rad) between 2 vectors v1 and v2 2071 ! dist = atan2 ( |v1 x v2| / (v1 . v2) ) 2072 ! (more accurate than acos(v1. v2) when dist close to 0 or Pi. 2073 ! 2011-08-25: replaced ACOS with ATAN2 2074 !======================================================================= 2075 real(kind=DP), intent(IN), dimension(1:) :: v1, v2 2076 real(kind=DP), intent(OUT) :: dist 2077 2078 real(kind=DP), dimension(1:3) :: v3 2079 real(kind=DP) :: sprod, vprod 2080 !======================================================================= 2081 2082 ! scalar product s = A. cos(theta) 2083 sprod = dot_product(v1, v2) 2084 ! vectorial product v = A. sin(theta) 2085 call vect_prod(v1, v2, v3) 2086 vprod = sqrt(dot_product(v3,v3)) 2087 ! theta = atan( |v|/s) in [0,Pi] 2088 dist = atan2( vprod, sprod) 2089 2090 return 2091 end subroutine angdist 2092 !======================================================================= 2093 subroutine vect_prod(v1, v2, v3) ! 2094 !======================================================================= 2095 ! returns in v3 the vectorial product of the 2 3-dimensional 2096 ! vectors v1 and v2 2097 !======================================================================= 2098 real(kind=DP), dimension(1:), INTENT(IN) :: v1, v2 2099 real(kind=DP), dimension(1:), INTENT(OUT) :: v3 2100 !======================================================================= 2101 2102 v3(1) = v1(2) * v2(3) - v1(3) * v2(2) 2103 v3(2) = v1(3) * v2(1) - v1(1) * v2(3) 2104 v3(3) = v1(1) * v2(2) - v1(2) * v2(1) 2105 2106 return 2107 end subroutine vect_prod 2108 2109 !**************************************************************************** 2110 2111 !======================================================================= 2112 function nside2ntemplates(nside) result(ntemplates) 2113 !======================================================================= 2114 ! returns the number of template pixels 2115 ! 2009-03-24: accepts Nside > 8192, and returns I8B result 2116 ! 2117 2118 integer(kind=i8b) :: ntemplates 2119 integer(kind=i4b), intent(IN) :: nside 2120 2121 ntemplates = 1_I8B + nside * (nside + 6_I8B) 2122 ntemplates = ntemplates / 4_I8B 2123 2124 return 2125 end function nside2ntemplates 2126 2127!****************************************************************************** 2128 function nside2npweights(nside) result(npweights) 2129!--------------------------------------- 2130! returns the number of full sky weights (in compressed form) for a given Nside. 2131! Given the symmetries of the Healpix lay out, the number of non-redundant weights 2132! is nf = ((3*Nside+1)*(Nside+1))/4 ~ Npix/16 2133! 2134! adapted from nside2npweights.pro 2135! 2018-05-18 2136!--------------------------------------- 2137 2138 integer(i8b) :: npweights 2139 integer(i4b), intent(in) :: nside 2140 integer(i8b), parameter :: one = 1_i8b 2141 integer(i8b), parameter :: four = 4_i8b 2142 2143 npweights = ((3*Nside + one)*(Nside + one)) / four 2144 2145 return 2146 end function nside2npweights 2147 2148 2149end module pix_tools 2150 2151 2152 2153