1; ----------------------------------------------------------------------------- 2; 3; Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, Anthony J. Banday 4; 5; 6; 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; ----------------------------------------------------------------------------- 28 29;=========================================================================== 30 31pro query_strip, nside, theta1, theta2, listpix, nlist, help=help, nested=nested, inclusive=inclusive 32;+ 33;======================================================================= 34; 35; query_strip, Nside, theta1, theta2, Listpix, Nlist, HELP=, NESTED=, INCLUSIVE= 36; -------------- 37; nside = resolution parameter (a power of 2) 38; theta1,theta2 = colatitude in Rad in [0,Pi] 39; if theta1<theta2 then returns pixels with 40; theta1<colat<theta2 41; if theta2<theta1 then returns 42; 0<= colatitude < theta2 or theta1 < colatitude < Pi 43; 44; list_pix = list of pixel lying in the strip 45; nlist = number of pixels in the list 46; nested (OPT), :0 by default, the output list is in RING scheme 47; if set to 1, the output list is in NESTED scheme 48; inclusive (OPT) , :0 by default, only the pixels whose center 49; lie in the strip are listed on output 50; if set to 1, all pixels overlapping the strip are output 51; help (OPT): prints this documentation header and exits 52; 53; 54; 55; v1.0, EH, IAP, 2008-03-28 : adapted from F90 code 56; 2012-01-14: systematically returns Listpix=[-1], Nlist=0 in case of problem 57;======================================================================= 58;- 59 60routine = 'query_strip' 61syntax = routine+', Nside, theta1, theta2, Listpix [, Nlist, HELP=, NESTED=, INCLUSIVE=]' 62nlist = 0 & listpix = [-1] 63 64if keyword_set(help) then begin 65 doc_library,routine 66 return 67endif 68 69if (n_params() lt 4 or n_params() gt 5) then begin 70 print,syntax 71 return 72endif 73 74 75npix = nside2npix(nside) 76if (npix lt 0) then begin 77 message,"Invalid Nside = "+string(nside) 78endif 79lnside = long(nside) 80 81if (theta1 lt 0.d0 or theta1 gt !PI or theta2 lt 0.d0 or theta2 gt !PI) then begin 82 message,/info," the colatitudes are in RADIAN " 83 message,/info," and should lie in [0,Pi] " 84 message,/info," current value = ", theta1, theta2 85 message," program abort " 86endif 87 88if (theta1 lt theta2) then begin 89 nstrip = 1 90 colrange = [ theta1, theta2 ] 91endif else begin 92 nstrip = 2 93 colrange = [ 0.0d0, theta2, theta1, !DPI ] 94endelse 95 96do_inclusive = keyword_set(inclusive) 97 98 99nlist=-1 100for is =0, nstrip-1 do begin 101 zu = cos(colrange[2*is ]) ; upper (=Northern) bound in z 102 zd = cos(colrange[2*is+1]) ; lower (=Southern) bound in z 103 irmin = ring_num(lnside, zu, shift=-do_inclusive) ; shift up north for inclusive 104 irmax = ring_num(lnside, zd, shift= do_inclusive) ; shift down south for inclusive 105 106 107 ; ------------- loop on ring number --------------------- 108 for iz = irmin, irmax do begin 109 110 ; z of ring being considered 111 zring = ring2z(lnside, iz) 112 113 if ((zring ge zd and zring le zu) or do_inclusive) then begin 114 ; ------- finds pixels in the ring --------- 115 listir=in_ring(lnside, iz, 0.d0, !DPI, nir, nested=nested) 116 117 if nlist le 0 then begin 118 listpix = listir 119 nlist = nir 120 endif else begin 121 listpix = [listpix,listir] 122 nlist = nlist + nir 123 endelse 124 125 endif 126 127 endfor 128 129endfor 130 131 132return 133end 134 135