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