Line data Source code
1 : /*
2 : * This file is part of Healpix_cxx.
3 : *
4 : * Healpix_cxx is free software; you can redistribute it and/or modify
5 : * it under the terms of the GNU General Public License as published by
6 : * the Free Software Foundation; either version 2 of the License, or
7 : * (at your option) any later version.
8 : *
9 : * Healpix_cxx is distributed in the hope that it will be useful,
10 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : * GNU General Public License for more details.
13 : *
14 : * You should have received a copy of the GNU General Public License
15 : * along with Healpix_cxx; if not, write to the Free Software
16 : * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17 : *
18 : * For more information about HEALPix, see http://healpix.sourceforge.net
19 : */
20 :
21 : /*
22 : * Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
23 : * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
24 : * (DLR).
25 : */
26 :
27 : /*
28 : * Copyright (C) 2003-2012 Max-Planck-Society
29 : * Author: Martin Reinecke
30 : */
31 :
32 : #include "healpix_base/healpix_base.h"
33 : #include "healpix_base/geom_utils.h"
34 : #include "healpix_base/lsconstants.h"
35 :
36 : namespace healpix{
37 :
38 : using namespace std;
39 :
40 : template<> const int T_Healpix_Base<int >::order_max=13;
41 : template<> const int T_Healpix_Base<int64>::order_max=29;
42 :
43 0 : template<typename I> int T_Healpix_Base<I>::nside2order (I nside)
44 : {
45 0 : planck_assert (nside>I(0), "invalid value for Nside");
46 0 : return ((nside)&(nside-1)) ? -1 : ilog2(nside);
47 : }
48 0 : template<typename I> I T_Healpix_Base<I>::npix2nside (I npix)
49 : {
50 0 : I res=isqrt(npix/I(12));
51 0 : planck_assert (npix==res*res*I(12), "invalid value for npix");
52 0 : return res;
53 : }
54 :
55 0 : template<typename I> I T_Healpix_Base<I>::ring_above (double z) const
56 : {
57 : double az=abs(z);
58 0 : if (az<=twothird) // equatorial region
59 0 : return I(nside_*(2-1.5*z));
60 0 : I iring = I(nside_*sqrt(3*(1-az)));
61 0 : return (z>0) ? iring : 4*nside_-iring-1;
62 : }
63 :
64 : namespace {
65 :
66 : /* Short note on the "zone":
67 : zone = 0: pixel lies completely outside the queried shape
68 : 1: pixel may overlap with the shape, pixel center is outside
69 : 2: pixel center is inside the shape, but maybe not the complete pixel
70 : 3: pixel lies completely inside the shape */
71 :
72 0 : template<typename I, typename I2> inline void check_pixel (int o, int order_,
73 : int omax, int zone, rangeset<I2> &pixset, I pix, vector<pair<I,int> > &stk,
74 : bool inclusive, int &stacktop)
75 : {
76 0 : if (zone==0) return;
77 :
78 0 : if (o<order_)
79 : {
80 0 : if (zone>=3)
81 : {
82 0 : int sdist=2*(order_-o); // the "bit-shift distance" between map orders
83 0 : pixset.append(pix<<sdist,(pix+1)<<sdist); // output all subpixels
84 : }
85 : else // (zone>=1)
86 0 : for (int i=0; i<4; ++i)
87 0 : stk.push_back(make_pair(4*pix+3-i,o+1)); // add children
88 : }
89 0 : else if (o>order_) // this implies that inclusive==true
90 : {
91 0 : if (zone>=2) // pixel center in shape
92 : {
93 0 : pixset.append(pix>>(2*(o-order_))); // output the parent pixel at order_
94 0 : stk.resize(stacktop); // unwind the stack
95 : }
96 : else // (zone>=1): pixel center in safety range
97 : {
98 0 : if (o<omax) // check sublevels
99 0 : for (int i=0; i<4; ++i) // add children in reverse order
100 0 : stk.push_back(make_pair(4*pix+3-i,o+1));
101 : else // at resolution limit
102 : {
103 0 : pixset.append(pix>>(2*(o-order_))); // output the parent pixel at order_
104 0 : stk.resize(stacktop); // unwind the stack
105 : }
106 : }
107 : }
108 : else // o==order_
109 : {
110 0 : if (zone>=2)
111 0 : pixset.append(pix);
112 0 : else if (inclusive) // and (zone>=1)
113 : {
114 0 : if (order_<omax) // check sublevels
115 : {
116 0 : stacktop=stk.size(); // remember current stack position
117 0 : for (int i=0; i<4; ++i) // add children in reverse order
118 0 : stk.push_back(make_pair(4*pix+3-i,o+1));
119 : }
120 : else // at resolution limit
121 0 : pixset.append(pix); // output the pixel
122 : }
123 : }
124 : }
125 :
126 0 : template<typename I> bool check_pixel_ring (const T_Healpix_Base<I> &b1,
127 : const T_Healpix_Base<I> &b2, I pix, I nr, I ipix1, int fct,
128 : double cz, double cphi, double cosrp2, I cpix)
129 : {
130 0 : if (pix>=nr) pix-=nr;
131 0 : if (pix<0) pix+=nr;
132 0 : pix+=ipix1;
133 0 : if (pix==cpix) return false; // disk center in pixel => overlap
134 : int px,py,pf;
135 0 : b1.pix2xyf(pix,px,py,pf);
136 0 : for (int i=0; i<fct-1; ++i) // go along the 4 edges
137 : {
138 0 : I ox=fct*px, oy=fct*py;
139 : double pz,pphi;
140 0 : b2.pix2zphi(b2.xyf2pix(ox+i,oy,pf),pz,pphi);
141 0 : if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
142 0 : return false;
143 0 : b2.pix2zphi(b2.xyf2pix(ox+fct-1,oy+i,pf),pz,pphi);
144 0 : if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
145 : return false;
146 0 : b2.pix2zphi(b2.xyf2pix(ox+fct-1-i,oy+fct-1,pf),pz,pphi);
147 0 : if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
148 : return false;
149 0 : b2.pix2zphi(b2.xyf2pix(ox,oy+fct-1-i,pf),pz,pphi);
150 0 : if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
151 : return false;
152 : }
153 : return true;
154 : }
155 :
156 : } // unnamed namespace
157 :
158 : template<typename I> template<typename I2>
159 0 : void T_Healpix_Base<I>::query_disc_internal
160 : (pointing ptg, double radius, int fact, rangeset<I2> &pixset) const
161 : {
162 0 : bool inclusive = (fact!=0);
163 : pixset.clear();
164 0 : ptg.normalize();
165 :
166 0 : if (scheme_==RING)
167 : {
168 : I fct=1;
169 0 : if (inclusive)
170 : {
171 0 : planck_assert (((I(1)<<order_max)/nside_)>=fact,
172 : "invalid oversampling factor");
173 : fct = fact;
174 : }
175 0 : T_Healpix_Base b2;
176 : double rsmall, rbig;
177 0 : if (fct>1)
178 : {
179 0 : b2.SetNside(fct*nside_,RING);
180 0 : rsmall = radius+b2.max_pixrad();
181 0 : rbig = radius+max_pixrad();
182 : }
183 : else
184 0 : rsmall = rbig = inclusive ? radius+max_pixrad() : radius;
185 :
186 0 : if (rsmall>=pi)
187 0 : { pixset.append(0,npix_); return; }
188 :
189 0 : rbig = min(pi,rbig);
190 :
191 0 : double cosrsmall = cos(rsmall);
192 0 : double cosrbig = cos(rbig);
193 :
194 0 : double z0 = cos(ptg.theta);
195 0 : double xa = 1./sqrt((1-z0)*(1+z0));
196 :
197 0 : I cpix=zphi2pix(z0,ptg.phi);
198 :
199 0 : double rlat1 = ptg.theta - rsmall;
200 0 : double zmax = cos(rlat1);
201 0 : I irmin = ring_above (zmax)+1;
202 :
203 0 : if ((rlat1<=0) && (irmin>1)) // north pole in the disk
204 : {
205 : I sp,rp; bool dummy;
206 0 : get_ring_info_small(irmin-1,sp,rp,dummy);
207 0 : pixset.append(0,sp+rp);
208 : }
209 :
210 0 : if ((fct>1) && (rlat1>0)) irmin=max(I(1),irmin-1);
211 :
212 0 : double rlat2 = ptg.theta + rsmall;
213 0 : double zmin = cos(rlat2);
214 0 : I irmax = ring_above (zmin);
215 :
216 0 : if ((fct>1) && (rlat2<pi)) irmax=min(4*nside_-1,irmax+1);
217 :
218 0 : for (I iz=irmin; iz<=irmax; ++iz)
219 : {
220 0 : double z=ring2z(iz);
221 0 : double x = (cosrbig-z*z0)*xa;
222 0 : double ysq = 1-z*z-x*x;
223 0 : double dphi = (ysq<=0) ? pi-1e-15 : atan2(sqrt(ysq),x);
224 : I nr, ipix1;
225 : bool shifted;
226 0 : get_ring_info_small(iz,ipix1,nr,shifted);
227 0 : double shift = shifted ? 0.5 : 0.;
228 :
229 0 : I ipix2 = ipix1 + nr - 1; // highest pixel number in the ring
230 :
231 0 : I ip_lo = ifloor<I>(nr*inv_twopi*(ptg.phi-dphi) - shift)+1;
232 0 : I ip_hi = ifloor<I>(nr*inv_twopi*(ptg.phi+dphi) - shift);
233 :
234 0 : if (fct>1)
235 : {
236 0 : while ((ip_lo<=ip_hi) && check_pixel_ring
237 0 : (*this,b2,ip_lo,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
238 0 : ++ip_lo;
239 0 : while ((ip_hi>ip_lo) && check_pixel_ring
240 0 : (*this,b2,ip_hi,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
241 0 : --ip_hi;
242 : }
243 :
244 0 : if (ip_lo<=ip_hi)
245 : {
246 0 : if (ip_hi>=nr)
247 0 : { ip_lo-=nr; ip_hi-=nr; }
248 0 : if (ip_lo<0)
249 : {
250 0 : pixset.append(ipix1,ipix1+ip_hi+1);
251 0 : pixset.append(ipix1+ip_lo+nr,ipix2+1);
252 : }
253 : else
254 0 : pixset.append(ipix1+ip_lo,ipix1+ip_hi+1);
255 : }
256 : }
257 0 : if ((rlat2>=pi) && (irmax+1<4*nside_)) // south pole in the disk
258 : {
259 : I sp,rp; bool dummy;
260 0 : get_ring_info_small(irmax+1,sp,rp,dummy);
261 0 : pixset.append(sp,npix_);
262 : }
263 : }
264 : else // scheme_==NEST
265 : {
266 0 : if (radius>=pi) // disk covers the whole sphere
267 0 : { pixset.append(0,npix_); return; }
268 :
269 : int oplus = 0;
270 0 : if (inclusive)
271 : {
272 0 : planck_assert ((I(1)<<(order_max-order_))>=fact,
273 : "invalid oversampling factor");
274 0 : planck_assert ((fact&(fact-1))==0,
275 : "oversampling factor must be a power of 2");
276 0 : oplus=ilog2(fact);
277 : }
278 0 : int omax=order_+oplus; // the order up to which we test
279 :
280 : vec3 vptg(ptg);
281 0 : arr<T_Healpix_Base<I> > base(omax+1);
282 : arr<double> crpdr(omax+1), crmdr(omax+1);
283 0 : double cosrad=cos(radius);
284 0 : for (int o=0; o<=omax; ++o) // prepare data at the required orders
285 : {
286 0 : base[o].Set(o,NEST);
287 0 : double dr=base[o].max_pixrad(); // safety distance
288 0 : crpdr[o] = (radius+dr>pi) ? -1. : cos(radius+dr);
289 0 : crmdr[o] = (radius-dr<0.) ? 1. : cos(radius-dr);
290 : }
291 : vector<pair<I,int> > stk; // stack for pixel numbers and their orders
292 0 : stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
293 0 : for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
294 0 : stk.push_back(make_pair(I(11-i),0));
295 :
296 0 : int stacktop=0; // a place to save a stack position
297 :
298 0 : while (!stk.empty()) // as long as there are pixels on the stack
299 : {
300 : // pop current pixel number and order from the stack
301 0 : I pix=stk.back().first;
302 0 : int o=stk.back().second;
303 : stk.pop_back();
304 :
305 : double z,phi;
306 : base[o].pix2zphi(pix,z,phi);
307 : // cosine of angular distance between pixel center and disk center
308 0 : double cangdist=cosdist_zphi(vptg.z,ptg.phi,z,phi);
309 :
310 0 : if (cangdist>crpdr[o])
311 : {
312 0 : int zone = (cangdist<cosrad) ? 1 : ((cangdist<=crmdr[o]) ? 2 : 3);
313 :
314 0 : check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
315 : stacktop);
316 : }
317 : }
318 : }
319 : }
320 :
321 0 : template<typename I> void T_Healpix_Base<I>::query_disc
322 : (pointing ptg, double radius, rangeset<I> &pixset) const
323 : {
324 0 : query_disc_internal (ptg, radius, 0, pixset);
325 0 : }
326 :
327 0 : template<typename I> void T_Healpix_Base<I>::query_disc_inclusive
328 : (pointing ptg, double radius, rangeset<I> &pixset, int fact) const
329 : {
330 0 : planck_assert(fact>0,"fact must be a positive integer");
331 0 : if ((sizeof(I)<8) && (((I(1)<<order_max)/nside_)<fact))
332 : {
333 0 : T_Healpix_Base<int64> base2(nside_,scheme_,SET_NSIDE);
334 0 : base2.query_disc_internal(ptg,radius,fact,pixset);
335 : return;
336 : }
337 0 : query_disc_internal (ptg, radius, fact, pixset);
338 : }
339 :
340 : template<typename I> template<typename I2>
341 0 : void T_Healpix_Base<I>::query_multidisc (const arr<vec3> &norm,
342 : const arr<double> &rad, int fact, rangeset<I2> &pixset) const
343 : {
344 0 : bool inclusive = (fact!=0);
345 : tsize nv=norm.size();
346 0 : planck_assert(nv==rad.size(),"inconsistent input arrays");
347 : pixset.clear();
348 :
349 0 : if (scheme_==RING)
350 : {
351 : I fct=1;
352 0 : if (inclusive)
353 : {
354 0 : planck_assert (((I(1)<<order_max)/nside_)>=fact,
355 : "invalid oversampling factor");
356 : fct = fact;
357 : }
358 0 : T_Healpix_Base b2;
359 : double rpsmall, rpbig;
360 0 : if (fct>1)
361 : {
362 0 : b2.SetNside(fct*nside_,RING);
363 0 : rpsmall = b2.max_pixrad();
364 0 : rpbig = max_pixrad();
365 : }
366 : else
367 0 : rpsmall = rpbig = inclusive ? max_pixrad() : 0;
368 :
369 0 : I irmin=1, irmax=4*nside_-1;
370 : vector<double> z0,xa,cosrsmall,cosrbig;
371 : vector<pointing> ptg;
372 : vector<I> cpix;
373 0 : for (tsize i=0; i<nv; ++i)
374 : {
375 0 : double rsmall=rad[i]+rpsmall;
376 0 : if (rsmall<pi)
377 : {
378 0 : double rbig=min(pi,rad[i]+rpbig);
379 : pointing pnt=pointing(norm[i]);
380 0 : cosrsmall.push_back(cos(rsmall));
381 0 : cosrbig.push_back(cos(rbig));
382 0 : double cth=cos(pnt.theta);
383 0 : z0.push_back(cth);
384 0 : if (fct>1) cpix.push_back(zphi2pix(cth,pnt.phi));
385 0 : xa.push_back(1./sqrt((1-cth)*(1+cth)));
386 0 : ptg.push_back(pnt);
387 :
388 0 : double rlat1 = pnt.theta - rsmall;
389 0 : double zmax = cos(rlat1);
390 0 : I irmin_t = (rlat1<=0) ? 1 : ring_above (zmax)+1;
391 :
392 0 : if ((fct>1) && (rlat1>0)) irmin_t=max(I(1),irmin_t-1);
393 :
394 0 : double rlat2 = pnt.theta + rsmall;
395 0 : double zmin = cos(rlat2);
396 0 : I irmax_t = (rlat2>=pi) ? 4*nside_-1 : ring_above (zmin);
397 :
398 0 : if ((fct>1) && (rlat2<pi)) irmax_t=min(4*nside_-1,irmax_t+1);
399 :
400 : if (irmax_t < irmax) irmax=irmax_t;
401 : if (irmin_t > irmin) irmin=irmin_t;
402 : }
403 : }
404 :
405 0 : for (I iz=irmin; iz<=irmax; ++iz)
406 : {
407 0 : double z=ring2z(iz);
408 : I ipix1,nr;
409 : bool shifted;
410 0 : get_ring_info_small(iz,ipix1,nr,shifted);
411 0 : double shift = shifted ? 0.5 : 0.;
412 : rangeset<I2> tr;
413 0 : tr.append(ipix1,ipix1+nr);
414 0 : for (tsize j=0; j<z0.size(); ++j)
415 : {
416 0 : double x = (cosrbig[j]-z*z0[j])*xa[j];
417 0 : double ysq = 1.-z*z-x*x;
418 0 : double dphi = (ysq<=0) ? pi-1e-15 : atan2(sqrt(ysq),x);
419 0 : I ip_lo = ifloor<I>(nr*inv_twopi*(ptg[j].phi-dphi) - shift)+1;
420 0 : I ip_hi = ifloor<I>(nr*inv_twopi*(ptg[j].phi+dphi) - shift);
421 0 : if (fct>1)
422 : {
423 0 : while ((ip_lo<=ip_hi) && check_pixel_ring
424 0 : (*this,b2,ip_lo,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
425 0 : ++ip_lo;
426 0 : while ((ip_hi>ip_lo) && check_pixel_ring
427 0 : (*this,b2,ip_hi,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
428 0 : --ip_hi;
429 : }
430 0 : if (ip_hi>=nr)
431 0 : { ip_lo-=nr; ip_hi-=nr;}
432 0 : if (ip_lo<0)
433 0 : tr.remove(ipix1+ip_hi+1,ipix1+ip_lo+nr);
434 : else
435 0 : tr.intersect(ipix1+ip_lo,ipix1+ip_hi+1);
436 : }
437 0 : pixset.append(tr);
438 : }
439 : }
440 : else // scheme_ == NEST
441 : {
442 : int oplus = 0;
443 0 : if (inclusive)
444 : {
445 0 : planck_assert ((I(1)<<(order_max-order_))>=fact,
446 : "invalid oversampling factor");
447 0 : planck_assert ((fact&(fact-1))==0,
448 : "oversampling factor must be a power of 2");
449 0 : oplus=ilog2(fact);
450 : }
451 0 : int omax=order_+oplus; // the order up to which we test
452 :
453 : // TODO: ignore all disks with radius>=pi
454 :
455 0 : arr<T_Healpix_Base<I> > base(omax+1);
456 : arr3<double> crlimit(omax+1,nv,3);
457 0 : for (int o=0; o<=omax; ++o) // prepare data at the required orders
458 : {
459 0 : base[o].Set(o,NEST);
460 0 : double dr=base[o].max_pixrad(); // safety distance
461 0 : for (tsize i=0; i<nv; ++i)
462 : {
463 0 : crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
464 0 : crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
465 0 : crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
466 : }
467 : }
468 :
469 : vector<pair<I,int> > stk; // stack for pixel numbers and their orders
470 0 : stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
471 0 : for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
472 0 : stk.push_back(make_pair(I(11-i),0));
473 :
474 0 : int stacktop=0; // a place to save a stack position
475 :
476 0 : while (!stk.empty()) // as long as there are pixels on the stack
477 : {
478 : // pop current pixel number and order from the stack
479 0 : I pix=stk.back().first;
480 0 : int o=stk.back().second;
481 : stk.pop_back();
482 :
483 0 : vec3 pv(base[o].pix2vec(pix));
484 :
485 : tsize zone=3;
486 0 : for (tsize i=0; i<nv; ++i)
487 : {
488 : double crad=dotprod(pv,norm[i]);
489 0 : for (tsize iz=0; iz<zone; ++iz)
490 0 : if (crad<crlimit(o,i,iz))
491 0 : if ((zone=iz)==0) goto bailout;
492 : }
493 :
494 0 : check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
495 : stacktop);
496 0 : bailout:;
497 : }
498 : }
499 0 : }
500 :
501 0 : template<typename I> void T_Healpix_Base<I>::query_multidisc_general
502 : (const arr<vec3> &norm, const arr<double> &rad, bool inclusive,
503 : const vector<int> &cmds, rangeset<I> &pixset) const
504 : {
505 : tsize nv=norm.size();
506 0 : planck_assert(nv==rad.size(),"inconsistent input arrays");
507 : pixset.clear();
508 :
509 0 : if (scheme_==RING)
510 : {
511 0 : planck_fail ("not yet implemented");
512 : }
513 : else // scheme_ == NEST
514 : {
515 0 : int oplus=inclusive ? 2 : 0;
516 0 : int omax=min(order_max,order_+oplus); // the order up to which we test
517 :
518 : // TODO: ignore all disks with radius>=pi
519 :
520 0 : arr<T_Healpix_Base<I> > base(omax+1);
521 : arr3<double> crlimit(omax+1,nv,3);
522 0 : for (int o=0; o<=omax; ++o) // prepare data at the required orders
523 : {
524 0 : base[o].Set(o,NEST);
525 0 : double dr=base[o].max_pixrad(); // safety distance
526 0 : for (tsize i=0; i<nv; ++i)
527 : {
528 0 : crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
529 0 : crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
530 0 : crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
531 : }
532 : }
533 :
534 : vector<pair<I,int> > stk; // stack for pixel numbers and their orders
535 0 : stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
536 0 : for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
537 0 : stk.push_back(make_pair(I(11-i),0));
538 :
539 0 : int stacktop=0; // a place to save a stack position
540 : arr<tsize> zone(nv);
541 :
542 0 : vector<tsize> zstk; zstk.reserve(cmds.size());
543 :
544 0 : while (!stk.empty()) // as long as there are pixels on the stack
545 : {
546 : // pop current pixel number and order from the stack
547 0 : I pix=stk.back().first;
548 0 : int o=stk.back().second;
549 : stk.pop_back();
550 :
551 0 : vec3 pv(base[o].pix2vec(pix));
552 :
553 0 : for (tsize i=0; i<nv; ++i)
554 : {
555 0 : zone[i]=3;
556 : double crad=dotprod(pv,norm[i]);
557 0 : for (tsize iz=0; iz<zone[i]; ++iz)
558 0 : if (crad<crlimit(o,i,iz))
559 0 : zone[i]=iz;
560 : }
561 :
562 0 : for (tsize i=0; i<cmds.size(); ++i)
563 : {
564 : tsize tmp;
565 0 : switch (cmds[i])
566 : {
567 : case -1: // union
568 0 : tmp=zstk.back(); zstk.pop_back();
569 0 : zstk.back() = max(zstk.back(),tmp);
570 0 : break;
571 : case -2: // intersection
572 0 : tmp=zstk.back(); zstk.pop_back();
573 0 : zstk.back() = min(zstk.back(),tmp);
574 0 : break;
575 : default: // add value
576 0 : zstk.push_back(zone[cmds[i]]);
577 : }
578 : }
579 0 : planck_assert(zstk.size()==1,"inconsistent commands");
580 0 : tsize zn=zstk[0]; zstk.pop_back();
581 :
582 0 : check_pixel (o, order_, omax, zn, pixset, pix, stk, inclusive,
583 : stacktop);
584 : }
585 : }
586 0 : }
587 :
588 : template<> inline int T_Healpix_Base<int>::spread_bits (int v) const
589 421 : { return utab[v&0xff] | (utab[(v>>8)&0xff]<<16); }
590 0 : template<> inline int64 T_Healpix_Base<int64>::spread_bits (int v) const
591 : {
592 0 : return int64(utab[ v &0xff]) | (int64(utab[(v>> 8)&0xff])<<16)
593 0 : | (int64(utab[(v>>16)&0xff])<<32) | (int64(utab[(v>>24)&0xff])<<48);
594 : }
595 :
596 : template<> inline int T_Healpix_Base<int>::compress_bits (int v) const
597 : {
598 0 : int raw = (v&0x5555) | ((v&0x55550000)>>15);
599 0 : return ctab[raw&0xff] | (ctab[raw>>8]<<4);
600 : }
601 842 : template<> inline int T_Healpix_Base<int64>::compress_bits (int64 v) const
602 : {
603 842 : int64 raw = v&0x5555555555555555ull;
604 842 : raw|=raw>>15;
605 842 : return ctab[ raw &0xff] | (ctab[(raw>> 8)&0xff]<< 4)
606 842 : | (ctab[(raw>>32)&0xff]<<16) | (ctab[(raw>>40)&0xff]<<20);
607 : }
608 :
609 421 : template<typename I> inline void T_Healpix_Base<I>::nest2xyf (I pix, int &ix,
610 : int &iy, int &face_num) const
611 : {
612 421 : face_num = pix>>(2*order_);
613 421 : pix &= (npface_-1);
614 421 : ix = compress_bits(pix);
615 421 : iy = compress_bits(pix>>1);
616 421 : }
617 :
618 421 : template<typename I> inline I T_Healpix_Base<I>::xyf2nest (int ix, int iy,
619 : int face_num) const
620 421 : { return (I(face_num)<<(2*order_)) + spread_bits(ix) + (spread_bits(iy)<<1); }
621 :
622 421 : template<typename I> void T_Healpix_Base<I>::ring2xyf (I pix, int &ix, int &iy,
623 : int &face_num) const
624 : {
625 : I iring, iphi, kshift, nr;
626 421 : I nl2 = 2*nside_;
627 :
628 421 : if (pix<ncap_) // North Polar cap
629 : {
630 0 : iring = (1+isqrt(1+2*pix))>>1; //counted from North pole
631 0 : iphi = (pix+1) - 2*iring*(iring-1);
632 : kshift = 0;
633 : nr = iring;
634 0 : face_num=(iphi-1)/nr;
635 : }
636 421 : else if (pix<(npix_-ncap_)) // Equatorial region
637 : {
638 421 : I ip = pix - ncap_;
639 421 : I tmp = (order_>=0) ? ip>>(order_+2) : ip/(4*nside_);
640 421 : iring = tmp+nside_;
641 421 : iphi = ip-tmp*4*nside_ + 1;
642 421 : kshift = (iring+nside_)&1;
643 : nr = nside_;
644 421 : I ire = iring-nside_+1,
645 421 : irm = nl2+2-ire;
646 421 : I ifm = iphi - ire/2 + nside_ -1,
647 421 : ifp = iphi - irm/2 + nside_ -1;
648 421 : if (order_>=0)
649 421 : { ifm >>= order_; ifp >>= order_; }
650 : else
651 0 : { ifm /= nside_; ifp /= nside_; }
652 421 : face_num = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8));
653 : }
654 : else // South Polar cap
655 : {
656 0 : I ip = npix_ - pix;
657 0 : iring = (1+isqrt(2*ip-1))>>1; //counted from South pole
658 0 : iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
659 : kshift = 0;
660 : nr = iring;
661 0 : iring = 2*nl2-iring;
662 0 : face_num = 8 + (iphi-1)/nr;
663 : }
664 :
665 421 : I irt = iring - (jrll[face_num]*nside_) + 1;
666 421 : I ipt = 2*iphi- jpll[face_num]*nr - kshift -1;
667 421 : if (ipt>=nl2) ipt-=8*nside_;
668 :
669 421 : ix = (ipt-irt) >>1;
670 421 : iy = (-ipt-irt) >>1;
671 421 : }
672 :
673 0 : template<typename I> I T_Healpix_Base<I>::xyf2ring (int ix, int iy,
674 : int face_num) const
675 : {
676 0 : I nl4 = 4*nside_;
677 0 : I jr = (jrll[face_num]*nside_) - ix - iy - 1;
678 :
679 : I nr, kshift, n_before;
680 :
681 : bool shifted;
682 0 : get_ring_info_small(jr,n_before,nr,shifted);
683 0 : nr>>=2;
684 0 : kshift=1-shifted;
685 0 : I jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
686 0 : planck_assert(jp<=4*nr,"must not happen");
687 0 : if (jp<1) jp+=nl4; // assumption: if this triggers, then nl4==4*nr
688 :
689 0 : return n_before + jp - 1;
690 : }
691 :
692 0 : template<typename I> T_Healpix_Base<I>::T_Healpix_Base ()
693 0 : : order_(-1), nside_(0), npface_(0), ncap_(0), npix_(0),
694 0 : fact1_(0), fact2_(0), scheme_(RING) {}
695 :
696 29 : template<typename I> void T_Healpix_Base<I>::Set (int order,
697 : Healpix_Ordering_Scheme scheme)
698 : {
699 29 : planck_assert ((order>=0)&&(order<=order_max), "bad order");
700 29 : order_ = order;
701 29 : nside_ = I(1)<<order;
702 29 : npface_ = nside_<<order_;
703 29 : ncap_ = (npface_-nside_)<<1;
704 29 : npix_ = 12*npface_;
705 29 : fact2_ = 4./npix_;
706 29 : fact1_ = (nside_<<1)*fact2_;
707 29 : scheme_ = scheme;
708 29 : }
709 0 : template<typename I> void T_Healpix_Base<I>::SetNside (I nside,
710 : Healpix_Ordering_Scheme scheme)
711 : {
712 0 : order_ = nside2order(nside);
713 0 : planck_assert ((scheme!=NEST) || (order_>=0),
714 : "SetNside: nside must be power of 2 for nested maps");
715 0 : nside_ = nside;
716 0 : npface_ = nside_*nside_;
717 0 : ncap_ = (npface_-nside_)<<1;
718 0 : npix_ = 12*npface_;
719 0 : fact2_ = 4./npix_;
720 0 : fact1_ = (nside_<<1)*fact2_;
721 0 : scheme_ = scheme;
722 0 : }
723 :
724 0 : template<typename I> double T_Healpix_Base<I>::ring2z (I ring) const
725 : {
726 0 : if (ring<nside_)
727 0 : return 1 - ring*ring*fact2_;
728 0 : if (ring <=3*nside_)
729 0 : return (2*nside_-ring)*fact1_;
730 0 : ring=4*nside_ - ring;
731 0 : return ring*ring*fact2_ - 1;
732 : }
733 :
734 0 : template<typename I> I T_Healpix_Base<I>::pix2ring (I pix) const
735 : {
736 0 : if (scheme_==RING)
737 : {
738 0 : if (pix<ncap_) // North Polar cap
739 0 : return (1+I(isqrt(1+2*pix)))>>1; // counted from North pole
740 0 : else if (pix<(npix_-ncap_)) // Equatorial region
741 0 : return (pix-ncap_)/(4*nside_) + nside_; // counted from North pole
742 : else // South Polar cap
743 0 : return 4*nside_-((1+I(isqrt(2*(npix_-pix)-1)))>>1);
744 : }
745 : else
746 : {
747 : int face_num, ix, iy;
748 0 : nest2xyf(pix,ix,iy,face_num);
749 0 : return (I(jrll[face_num])<<order_) - ix - iy - 1;
750 : }
751 : }
752 :
753 0 : template<typename I> I T_Healpix_Base<I>::nest2ring (I pix) const
754 : {
755 0 : planck_assert(order_>=0, "hierarchical map required");
756 : int ix, iy, face_num;
757 0 : nest2xyf (pix, ix, iy, face_num);
758 0 : return xyf2ring (ix, iy, face_num);
759 : }
760 :
761 421 : template<typename I> I T_Healpix_Base<I>::ring2nest (I pix) const
762 : {
763 421 : planck_assert(order_>=0, "hierarchical map required");
764 : int ix, iy, face_num;
765 421 : ring2xyf (pix, ix, iy, face_num);
766 421 : return xyf2nest (ix, iy, face_num);
767 : }
768 :
769 0 : template<typename I> inline I T_Healpix_Base<I>::nest_peano_helper
770 : (I pix, int dir) const
771 : {
772 0 : int face = pix>>(2*order_);
773 : I result = 0;
774 0 : uint8 path = peano_face2path[dir][face];
775 :
776 0 : for (int shift=2*order_-2; shift>=0; shift-=2)
777 : {
778 0 : uint8 spix = uint8((pix>>shift) & 0x3);
779 0 : result <<= 2;
780 0 : result |= peano_subpix[dir][path][spix];
781 0 : path=peano_subpath[dir][path][spix];
782 : }
783 :
784 0 : return result + (I(peano_face2face[dir][face])<<(2*order_));
785 : }
786 :
787 0 : template<typename I> I T_Healpix_Base<I>::nest2peano (I pix) const
788 0 : { return nest_peano_helper(pix,0); }
789 :
790 0 : template<typename I> I T_Healpix_Base<I>::peano2nest (I pix) const
791 0 : { return nest_peano_helper(pix,1); }
792 :
793 24585 : template<typename I> I T_Healpix_Base<I>::loc2pix (double z, double phi,
794 : double sth, bool have_sth) const
795 : {
796 : double za = abs(z);
797 24585 : double tt = fmodulo(phi*inv_halfpi,4.0); // in [0,4)
798 :
799 24585 : if (scheme_==RING)
800 : {
801 24585 : if (za<=twothird) // Equatorial region
802 : {
803 16535 : I nl4 = 4*nside_;
804 16535 : double temp1 = nside_*(0.5+tt);
805 16535 : double temp2 = nside_*z*0.75;
806 16535 : I jp = I(temp1-temp2); // index of ascending edge line
807 16535 : I jm = I(temp1+temp2); // index of descending edge line
808 :
809 : // ring number counted from z=2/3
810 16535 : I ir = nside_ + 1 + jp - jm; // in {1,2n+1}
811 16535 : I kshift = 1-(ir&1); // kshift=1 if ir even, 0 otherwise
812 :
813 16535 : I t1 = jp+jm-nside_+kshift+1+nl4+nl4;
814 16535 : I ip = (order_>0) ?
815 16535 : (t1>>1)&(nl4-1) : ((t1>>1)%nl4); // in {0,4n-1}
816 :
817 16535 : return ncap_ + (ir-1)*nl4 + ip;
818 : }
819 : else // North & South polar caps
820 : {
821 8050 : double tp = tt-I(tt);
822 8050 : double tmp = ((za<0.99)||(!have_sth)) ?
823 7809 : nside_*sqrt(3*(1-za)) :
824 241 : nside_*sth/sqrt((1.+za)/3.);
825 :
826 8050 : I jp = I(tp*tmp); // increasing edge line index
827 8050 : I jm = I((1.0-tp)*tmp); // decreasing edge line index
828 :
829 8050 : I ir = jp+jm+1; // ring number counted from the closest pole
830 8050 : I ip = I(tt*ir); // in {0,4*ir-1}
831 8050 : planck_assert((ip>=0)&&(ip<4*ir),"must not happen");
832 : //ip %= 4*ir;
833 :
834 8050 : return (z>0) ? 2*ir*(ir-1) + ip : npix_ - 2*ir*(ir+1) + ip;
835 : }
836 : }
837 : else // scheme_ == NEST
838 : {
839 0 : if (za<=twothird) // Equatorial region
840 : {
841 0 : double temp1 = nside_*(0.5+tt);
842 0 : double temp2 = nside_*(z*0.75);
843 0 : I jp = I(temp1-temp2); // index of ascending edge line
844 0 : I jm = I(temp1+temp2); // index of descending edge line
845 0 : I ifp = jp >> order_; // in {0,4}
846 0 : I ifm = jm >> order_;
847 0 : int face_num = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8));
848 :
849 0 : int ix = jm & (nside_-1),
850 0 : iy = nside_ - (jp & (nside_-1)) - 1;
851 0 : return xyf2nest(ix,iy,face_num);
852 : }
853 : else // polar region, za > 2/3
854 : {
855 0 : int ntt = min(3,int(tt));
856 0 : double tp = tt-ntt;
857 0 : double tmp = ((za<0.99)||(!have_sth)) ?
858 0 : nside_*sqrt(3*(1-za)) :
859 0 : nside_*sth/sqrt((1.+za)/3.);
860 :
861 0 : I jp = I(tp*tmp); // increasing edge line index
862 0 : I jm = I((1.0-tp)*tmp); // decreasing edge line index
863 0 : jp=min(jp,nside_-1); // for points too close to the boundary
864 0 : jm=min(jm,nside_-1);
865 0 : return (z>=0) ?
866 0 : xyf2nest(nside_-jm -1,nside_-jp-1,ntt) : xyf2nest(jp,jm,ntt+8);
867 : }
868 : }
869 : }
870 :
871 135593 : template<typename I> void T_Healpix_Base<I>::pix2loc (I pix, double &z,
872 : double &phi, double &sth, bool &have_sth) const
873 : {
874 135593 : have_sth=false;
875 135593 : if (scheme_==RING)
876 : {
877 135172 : if (pix<ncap_) // North Polar cap
878 : {
879 22080 : I iring = (1+I(isqrt(1+2*pix)))>>1; // counted from North pole
880 22080 : I iphi = (pix+1) - 2*iring*(iring-1);
881 :
882 22080 : double tmp=(iring*iring)*fact2_;
883 22080 : z = 1.0 - tmp;
884 22080 : if (z>0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=true; }
885 22080 : phi = (iphi-0.5) * halfpi/iring;
886 : }
887 113092 : else if (pix<(npix_-ncap_)) // Equatorial region
888 : {
889 91010 : I nl4 = 4*nside_;
890 91010 : I ip = pix - ncap_;
891 91010 : I tmp = (order_>=0) ? ip>>(order_+2) : ip/nl4;
892 91010 : I iring = tmp + nside_,
893 91010 : iphi = ip-nl4*tmp+1;
894 : // 1 if iring+nside is odd, 1/2 otherwise
895 91010 : double fodd = ((iring+nside_)&1) ? 1 : 0.5;
896 :
897 91010 : z = (2*nside_-iring)*fact1_;
898 91010 : phi = (iphi-fodd) * pi*0.75*fact1_;
899 : }
900 : else // South Polar cap
901 : {
902 22082 : I ip = npix_ - pix;
903 22082 : I iring = (1+I(isqrt(2*ip-1)))>>1; // counted from South pole
904 22082 : I iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
905 :
906 22082 : double tmp=(iring*iring)*fact2_;
907 22082 : z = tmp - 1.0;
908 22082 : if (z<-0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=true; }
909 22082 : phi = (iphi-0.5) * halfpi/iring;
910 : }
911 : }
912 : else
913 : {
914 : int face_num, ix, iy;
915 421 : nest2xyf(pix,ix,iy,face_num);
916 :
917 421 : I jr = (I(jrll[face_num])<<order_) - ix - iy - 1;
918 :
919 : I nr;
920 421 : if (jr<nside_)
921 : {
922 : nr = jr;
923 0 : double tmp=(nr*nr)*fact2_;
924 0 : z = 1 - tmp;
925 0 : if (z>0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=true; }
926 : }
927 421 : else if (jr > 3*nside_)
928 : {
929 0 : nr = nside_*4-jr;
930 0 : double tmp=(nr*nr)*fact2_;
931 0 : z = tmp - 1;
932 0 : if (z<-0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=true; }
933 : }
934 : else
935 : {
936 : nr = nside_;
937 421 : z = (2*nside_-jr)*fact1_;
938 : }
939 :
940 421 : I tmp=I(jpll[face_num])*nr+ix-iy;
941 421 : planck_assert(tmp<8*nr,"must not happen");
942 421 : if (tmp<0) tmp+=8*nr;
943 421 : phi = (nr==nside_) ? 0.75*halfpi*tmp*fact1_ :
944 0 : (0.5*halfpi*tmp)/nr;
945 : }
946 135593 : }
947 :
948 : template<typename I> template<typename I2>
949 0 : void T_Healpix_Base<I>::query_polygon_internal
950 : (const vector<pointing> &vertex, int fact, rangeset<I2> &pixset) const
951 : {
952 : bool inclusive = (fact!=0);
953 : tsize nv=vertex.size();
954 0 : tsize ncirc = inclusive ? nv+1 : nv;
955 0 : planck_assert(nv>=3,"not enough vertices in polygon");
956 : arr<vec3> vv(nv);
957 0 : for (tsize i=0; i<nv; ++i)
958 0 : vv[i]=vertex[i].to_vec3();
959 : arr<vec3> normal(ncirc);
960 : int flip=0;
961 0 : for (tsize i=0; i<nv; ++i)
962 : {
963 0 : normal[i]=crossprod(vv[i],vv[(i+1)%nv]);
964 0 : double hnd=dotprod(normal[i],vv[(i+2)%nv]);
965 0 : planck_assert(abs(hnd)>1e-10,"degenerate corner");
966 0 : if (i==0)
967 0 : flip = (hnd<0.) ? -1 : 1;
968 : else
969 0 : planck_assert(flip*hnd>0,"polygon is not convex");
970 0 : normal[i]*=flip/normal[i].Length();
971 : }
972 : arr<double> rad(ncirc,halfpi);
973 0 : if (inclusive)
974 : {
975 : double cosrad;
976 0 : find_enclosing_circle (vv, normal[nv], cosrad);
977 0 : rad[nv]=acos(cosrad);
978 : }
979 0 : query_multidisc(normal,rad,fact,pixset);
980 0 : }
981 :
982 0 : template<typename I> void T_Healpix_Base<I>::query_polygon
983 : (const vector<pointing> &vertex, rangeset<I> &pixset) const
984 : {
985 0 : query_polygon_internal(vertex, 0, pixset);
986 0 : }
987 :
988 0 : template<typename I> void T_Healpix_Base<I>::query_polygon_inclusive
989 : (const vector<pointing> &vertex, rangeset<I> &pixset, int fact) const
990 : {
991 0 : planck_assert(fact>0,"fact must be a positive integer");
992 0 : if ((sizeof(I)<8) && (((I(1)<<order_max)/nside_)<fact))
993 : {
994 0 : T_Healpix_Base<int64> base2(nside_,scheme_,SET_NSIDE);
995 0 : base2.query_polygon_internal(vertex,fact,pixset);
996 : return;
997 : }
998 0 : query_polygon_internal(vertex, fact, pixset);
999 : }
1000 :
1001 0 : template<typename I> void T_Healpix_Base<I>::query_strip_internal
1002 : (double theta1, double theta2, bool inclusive, rangeset<I> &pixset) const
1003 : {
1004 0 : if (scheme_==RING)
1005 : {
1006 0 : I ring1 = max(I(1),1+ring_above(cos(theta1))),
1007 0 : ring2 = min(4*nside_-1,ring_above(cos(theta2)));
1008 0 : if (inclusive)
1009 : {
1010 0 : ring1 = max(I(1),ring1-1);
1011 0 : ring2 = min(4*nside_-1,ring2+1);
1012 : }
1013 :
1014 : I sp1,rp1,sp2,rp2;
1015 : bool dummy;
1016 0 : get_ring_info_small(ring1,sp1,rp1,dummy);
1017 0 : get_ring_info_small(ring2,sp2,rp2,dummy);
1018 0 : I pix1 = sp1,
1019 0 : pix2 = sp2+rp2;
1020 0 : if (pix1<=pix2) pixset.append(pix1,pix2);
1021 : }
1022 : else
1023 0 : planck_fail("query_strip not yet implemented for NESTED");
1024 0 : }
1025 :
1026 0 : template<typename I> void T_Healpix_Base<I>::query_strip (double theta1,
1027 : double theta2, bool inclusive, rangeset<I> &pixset) const
1028 : {
1029 : pixset.clear();
1030 :
1031 0 : if (theta1<theta2)
1032 0 : query_strip_internal(theta1,theta2,inclusive,pixset);
1033 : else
1034 : {
1035 0 : query_strip_internal(0.,theta2,inclusive,pixset);
1036 : rangeset<I> ps2;
1037 0 : query_strip_internal(theta1,pi,inclusive,ps2);
1038 0 : pixset.append(ps2);
1039 : }
1040 0 : }
1041 :
1042 0 : template<typename I> inline void T_Healpix_Base<I>::get_ring_info_small
1043 : (I ring, I &startpix, I &ringpix, bool &shifted) const
1044 : {
1045 0 : if (ring < nside_)
1046 : {
1047 0 : shifted = true;
1048 0 : ringpix = 4*ring;
1049 0 : startpix = 2*ring*(ring-1);
1050 : }
1051 0 : else if (ring < 3*nside_)
1052 : {
1053 0 : shifted = ((ring-nside_) & 1) == 0;
1054 0 : ringpix = 4*nside_;
1055 0 : startpix = ncap_ + (ring-nside_)*ringpix;
1056 : }
1057 : else
1058 : {
1059 0 : shifted = true;
1060 0 : I nr= 4*nside_-ring;
1061 0 : ringpix = 4*nr;
1062 0 : startpix = npix_-2*nr*(nr+1);
1063 : }
1064 0 : }
1065 :
1066 0 : template<typename I> void T_Healpix_Base<I>::get_ring_info (I ring, I &startpix,
1067 : I &ringpix, double &costheta, double &sintheta, bool &shifted) const
1068 : {
1069 0 : I northring = (ring>2*nside_) ? 4*nside_-ring : ring;
1070 0 : if (northring < nside_)
1071 : {
1072 0 : double tmp = northring*northring*fact2_;
1073 0 : costheta = 1 - tmp;
1074 0 : sintheta = sqrt(tmp*(2-tmp));
1075 0 : ringpix = 4*northring;
1076 0 : shifted = true;
1077 0 : startpix = 2*northring*(northring-1);
1078 : }
1079 : else
1080 : {
1081 0 : costheta = (2*nside_-northring)*fact1_;
1082 0 : sintheta = sqrt((1+costheta)*(1-costheta));
1083 0 : ringpix = 4*nside_;
1084 0 : shifted = ((northring-nside_) & 1) == 0;
1085 0 : startpix = ncap_ + (northring-nside_)*ringpix;
1086 : }
1087 0 : if (northring != ring) // southern hemisphere
1088 : {
1089 0 : costheta = -costheta;
1090 0 : startpix = npix_ - startpix - ringpix;
1091 : }
1092 0 : }
1093 0 : template<typename I> void T_Healpix_Base<I>::get_ring_info2 (I ring,
1094 : I &startpix, I &ringpix, double &theta, bool &shifted) const
1095 : {
1096 0 : I northring = (ring>2*nside_) ? 4*nside_-ring : ring;
1097 0 : if (northring < nside_)
1098 : {
1099 0 : double tmp = northring*northring*fact2_;
1100 0 : double costheta = 1 - tmp;
1101 0 : double sintheta = sqrt(tmp*(2-tmp));
1102 0 : theta = atan2(sintheta,costheta);
1103 0 : ringpix = 4*northring;
1104 0 : shifted = true;
1105 0 : startpix = 2*northring*(northring-1);
1106 : }
1107 : else
1108 : {
1109 0 : theta = acos((2*nside_-northring)*fact1_);
1110 0 : ringpix = 4*nside_;
1111 0 : shifted = ((northring-nside_) & 1) == 0;
1112 0 : startpix = ncap_ + (northring-nside_)*ringpix;
1113 : }
1114 0 : if (northring != ring) // southern hemisphere
1115 : {
1116 0 : theta = pi-theta;
1117 0 : startpix = npix_ - startpix - ringpix;
1118 : }
1119 0 : }
1120 :
1121 0 : template<typename I> void T_Healpix_Base<I>::neighbors (I pix,
1122 : fix_arr<I,8> &result) const
1123 : {
1124 : int ix, iy, face_num;
1125 0 : (scheme_==RING) ?
1126 0 : ring2xyf(pix,ix,iy,face_num) : nest2xyf(pix,ix,iy,face_num);
1127 :
1128 0 : const I nsm1 = nside_-1;
1129 0 : if ((ix>0)&&(ix<nsm1)&&(iy>0)&&(iy<nsm1))
1130 : {
1131 0 : if (scheme_==RING)
1132 0 : for (int m=0; m<8; ++m)
1133 0 : result[m] = xyf2ring(ix+nb_xoffset[m],iy+nb_yoffset[m],face_num);
1134 : else
1135 : {
1136 0 : I fpix = I(face_num)<<(2*order_),
1137 0 : px0=spread_bits(ix ), py0=spread_bits(iy )<<1,
1138 0 : pxp=spread_bits(ix+1), pyp=spread_bits(iy+1)<<1,
1139 0 : pxm=spread_bits(ix-1), pym=spread_bits(iy-1)<<1;
1140 :
1141 0 : result[0] = fpix+pxm+py0; result[1] = fpix+pxm+pyp;
1142 0 : result[2] = fpix+px0+pyp; result[3] = fpix+pxp+pyp;
1143 0 : result[4] = fpix+pxp+py0; result[5] = fpix+pxp+pym;
1144 0 : result[6] = fpix+px0+pym; result[7] = fpix+pxm+pym;
1145 : }
1146 : }
1147 : else
1148 : {
1149 0 : for (int i=0; i<8; ++i)
1150 : {
1151 0 : int x=ix+nb_xoffset[i], y=iy+nb_yoffset[i];
1152 : int nbnum=4;
1153 0 : if (x<0)
1154 0 : { x+=nside_; nbnum-=1; }
1155 0 : else if (x>=nside_)
1156 0 : { x-=nside_; nbnum+=1; }
1157 0 : if (y<0)
1158 0 : { y+=nside_; nbnum-=3; }
1159 0 : else if (y>=nside_)
1160 0 : { y-=nside_; nbnum+=3; }
1161 :
1162 0 : int f = nb_facearray[nbnum][face_num];
1163 0 : if (f>=0)
1164 : {
1165 0 : int bits = nb_swaparray[nbnum][face_num>>2];
1166 0 : if (bits&1) x=nside_-x-1;
1167 0 : if (bits&2) y=nside_-y-1;
1168 0 : if (bits&4) std::swap(x,y);
1169 0 : result[i] = (scheme_==RING) ? xyf2ring(x,y,f) : xyf2nest(x,y,f);
1170 : }
1171 : else
1172 0 : result[i] = -1;
1173 : }
1174 : }
1175 0 : }
1176 :
1177 0 : template<typename I> void T_Healpix_Base<I>::get_interpol (const pointing &ptg,
1178 : fix_arr<I,4> &pix, fix_arr<double,4> &wgt) const
1179 : {
1180 0 : planck_assert((ptg.theta>=0)&&(ptg.theta<=pi),"invalid theta value");
1181 0 : double z = cos (ptg.theta);
1182 0 : I ir1 = ring_above(z);
1183 0 : I ir2 = ir1+1;
1184 : double theta1, theta2, w1, tmp, dphi;
1185 : I sp,nr;
1186 : bool shift;
1187 : I i1,i2;
1188 0 : if (ir1>0)
1189 : {
1190 0 : get_ring_info2 (ir1, sp, nr, theta1, shift);
1191 0 : dphi = twopi/nr;
1192 0 : tmp = (ptg.phi/dphi - .5*shift);
1193 0 : i1 = (tmp<0) ? I(tmp)-1 : I(tmp);
1194 0 : w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
1195 0 : i2 = i1+1;
1196 0 : if (i1<0) i1 +=nr;
1197 0 : if (i2>=nr) i2 -=nr;
1198 0 : pix[0] = sp+i1; pix[1] = sp+i2;
1199 0 : wgt[0] = 1-w1; wgt[1] = w1;
1200 : }
1201 0 : if (ir2<(4*nside_))
1202 : {
1203 0 : get_ring_info2 (ir2, sp, nr, theta2, shift);
1204 0 : dphi = twopi/nr;
1205 0 : tmp = (ptg.phi/dphi - .5*shift);
1206 0 : i1 = (tmp<0) ? I(tmp)-1 : I(tmp);
1207 0 : w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
1208 0 : i2 = i1+1;
1209 0 : if (i1<0) i1 +=nr;
1210 0 : if (i2>=nr) i2 -=nr;
1211 0 : pix[2] = sp+i1; pix[3] = sp+i2;
1212 0 : wgt[2] = 1-w1; wgt[3] = w1;
1213 : }
1214 :
1215 0 : if (ir1==0)
1216 : {
1217 0 : double wtheta = ptg.theta/theta2;
1218 0 : wgt[2] *= wtheta; wgt[3] *= wtheta;
1219 0 : double fac = (1-wtheta)*0.25;
1220 0 : wgt[0] = fac; wgt[1] = fac; wgt[2] += fac; wgt[3] +=fac;
1221 0 : pix[0] = (pix[2]+2)&3;
1222 0 : pix[1] = (pix[3]+2)&3;
1223 : }
1224 0 : else if (ir2==4*nside_)
1225 : {
1226 0 : double wtheta = (ptg.theta-theta1)/(pi-theta1);
1227 0 : wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
1228 0 : double fac = wtheta*0.25;
1229 0 : wgt[0] += fac; wgt[1] += fac; wgt[2] = fac; wgt[3] =fac;
1230 0 : pix[2] = ((pix[0]+2)&3)+npix_-4;
1231 0 : pix[3] = ((pix[1]+2)&3)+npix_-4;
1232 : }
1233 : else
1234 : {
1235 0 : double wtheta = (ptg.theta-theta1)/(theta2-theta1);
1236 0 : wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
1237 0 : wgt[2] *= wtheta; wgt[3] *= wtheta;
1238 : }
1239 :
1240 0 : if (scheme_==NEST)
1241 0 : for (tsize m=0; m<pix.size(); ++m)
1242 0 : pix[m] = ring2nest(pix[m]);
1243 0 : }
1244 :
1245 0 : template<typename I> void T_Healpix_Base<I>::swap (T_Healpix_Base &other)
1246 : {
1247 : std::swap(order_,other.order_);
1248 : std::swap(nside_,other.nside_);
1249 : std::swap(npface_,other.npface_);
1250 : std::swap(ncap_,other.ncap_);
1251 : std::swap(npix_,other.npix_);
1252 : std::swap(fact1_,other.fact1_);
1253 : std::swap(fact2_,other.fact2_);
1254 : std::swap(scheme_,other.scheme_);
1255 0 : }
1256 :
1257 0 : template<typename I> double T_Healpix_Base<I>::max_pixrad() const
1258 : {
1259 : vec3 va,vb;
1260 0 : va.set_z_phi (2./3., pi/(4*nside_));
1261 0 : double t1 = 1.-1./nside_;
1262 0 : t1*=t1;
1263 0 : vb.set_z_phi (1-t1/3, 0);
1264 0 : return v_angle(va,vb);
1265 : }
1266 :
1267 0 : template<typename I> double T_Healpix_Base<I>::max_pixrad(I ring) const
1268 : {
1269 0 : if (ring>=2*nside_) ring=4*nside_-ring;
1270 0 : double z=ring2z(ring), z_up=(ring>1) ? ring2z(ring-1) : 1.;
1271 : vec3 mypos, uppos;
1272 0 : uppos.set_z_phi(z_up,0);
1273 0 : if (ring<=nside_)
1274 : {
1275 0 : mypos.set_z_phi(z,pi/(4*ring));
1276 0 : return v_angle(mypos,uppos);
1277 : }
1278 0 : mypos.set_z_phi(z,0);
1279 0 : double vdist=v_angle(mypos,uppos);
1280 0 : double hdist=sqrt(1.-z*z)*pi/(4*nside_);
1281 0 : return max(hdist,vdist);
1282 : }
1283 :
1284 0 : template<typename I> void T_Healpix_Base<I>::xyf2loc (double x, double y,
1285 : int face, double &z, double &phi, double &sth, bool &have_sth) const
1286 : {
1287 0 : have_sth = false;
1288 0 : double jr = jrll[face] - x - y;
1289 : double nr;
1290 0 : if (jr<1)
1291 : {
1292 : nr = jr;
1293 0 : double tmp = nr*nr/3.;
1294 0 : z = 1 - tmp;
1295 0 : if (z > 0.99)
1296 : {
1297 0 : sth = std::sqrt(tmp*(2.0-tmp));
1298 0 : have_sth = true;
1299 : }
1300 : }
1301 0 : else if (jr>3)
1302 : {
1303 0 : nr = 4-jr;
1304 0 : double tmp = nr*nr/3.;
1305 0 : z = tmp - 1;
1306 0 : if (z<-0.99)
1307 : {
1308 0 : sth = std::sqrt(tmp*(2.-tmp));
1309 0 : have_sth = true;
1310 : }
1311 : }
1312 : else
1313 : {
1314 : nr = 1;
1315 0 : z = (2-jr)*2./3.;
1316 : }
1317 :
1318 0 : double tmp=jpll[face]*nr+x-y;
1319 0 : if (tmp<0) tmp+=8;
1320 0 : if (tmp>=8) tmp-=8;
1321 0 : phi = (nr<1e-15) ? 0 : (0.5*halfpi*tmp)/nr;
1322 0 : }
1323 :
1324 : namespace {
1325 :
1326 0 : vec3 locToVec3 (double z, double phi, double sth, bool have_sth)
1327 : {
1328 0 : if (have_sth)
1329 0 : return vec3(sth*cos(phi),sth*sin(phi),z);
1330 : else
1331 : {
1332 : vec3 res;
1333 0 : res.set_z_phi (z, phi);
1334 0 : return res;
1335 : }
1336 : }
1337 :
1338 : } // unnamed namespace
1339 :
1340 0 : template<typename I> void T_Healpix_Base<I>::boundaries(I pix, tsize step,
1341 : vector<vec3> &out) const
1342 : {
1343 0 : out.resize(4*step);
1344 : int ix, iy, face;
1345 0 : pix2xyf(pix, ix, iy, face);
1346 0 : double dc = 0.5 / nside_;
1347 0 : double xc = (ix + 0.5)/nside_, yc = (iy + 0.5)/nside_;
1348 0 : double d = 1.0/(step*nside_);
1349 0 : for (tsize i=0; i<step; ++i)
1350 : {
1351 : double z, phi, sth;
1352 : bool have_sth;
1353 0 : xyf2loc(xc+dc-i*d, yc+dc, face, z, phi, sth, have_sth);
1354 0 : out[i] = locToVec3(z, phi, sth, have_sth);
1355 0 : xyf2loc(xc-dc, yc+dc-i*d, face, z, phi, sth, have_sth);
1356 0 : out[i+step] = locToVec3(z, phi, sth, have_sth);
1357 0 : xyf2loc(xc-dc+i*d, yc-dc, face, z, phi, sth, have_sth);
1358 0 : out[i+2*step] = locToVec3(z, phi, sth, have_sth);
1359 0 : xyf2loc(xc+dc, yc-dc+i*d, face, z, phi, sth, have_sth);
1360 0 : out[i+3*step] = locToVec3(z, phi, sth, have_sth);
1361 : }
1362 0 : }
1363 :
1364 0 : template<typename I> arr<int> T_Healpix_Base<I>::swap_cycles() const
1365 : {
1366 0 : planck_assert(order_>=0, "need hierarchical map");
1367 0 : planck_assert(order_<=13, "map too large");
1368 0 : arr<int> result(swap_clen[order_]);
1369 : tsize ofs=0;
1370 0 : for (int m=0; m<order_;++m) ofs+=swap_clen[m];
1371 0 : for (tsize m=0; m<result.size();++m) result[m]=swap_cycle[m+ofs];
1372 0 : return result;
1373 : }
1374 :
1375 : template class T_Healpix_Base<int>;
1376 : template class T_Healpix_Base<int64>;
1377 : } // namespace healpix
1378 :
|
