source: Sophya/trunk/SophyaLib/SkyMap/HEALPixUtils.cc@ 1195

Last change on this file since 1195 was 1195, checked in by ansari, 25 years ago

fichier d'utilitaire HEALPix
File size: 19.5 KB
Line 
1// utilitaires de pixelisation HEALPix
2#include "HEALPixUtils.h"
3#include <iostream.h>
4#include <math.h>
5//#include <complex>
6#include "tvector.h"
7#include "smathconst.h"
8extern "C"
9{
10#include <stdio.h>
11#include <stdlib.h>
12#include <unistd.h>
13}
14
15namespace SOPHYA
16{
17
18
19
20//*******************************************************************
21// Class PIXELS_XY
22// Construction des tableaux necessaires a la traduction des indices RING en
23// indices NESTED (ou l'inverse)
24//*******************************************************************
25
26PIXELS_XY::PIXELS_XY()
27{
28 pix2x_.ReSize(1024);
29 pix2x_.Reset();
30 pix2y_.ReSize(1024);
31 pix2y_.Reset();
32 x2pix_.ReSize(128);
33 x2pix_.Reset();
34 y2pix_.ReSize(128);
35 y2pix_.Reset();
36 mk_pix2xy();
37 mk_xy2pix();
38}
39
40PIXELS_XY& PIXELS_XY::instance()
41{
42 static PIXELS_XY single;
43 return (single);
44}
45
46void PIXELS_XY::mk_pix2xy()
47{
48 /*
49 ==================================================
50 subroutine mk_pix2xy
51 ==================================================
52 c constructs the array giving x and y in the face from pixel number
53 c for the nested (quad-cube like) ordering of pixels
54 c
55 c the bits corresponding to x and y are interleaved in the pixel number
56 c one breaks up the pixel number by even and odd bits
57 ==================================================
58 */
59 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
60 // (16/12/98)
61
62 int kpix, jpix, IX, IY, IP, ID;
63
64 for(kpix = 0; kpix < 1024; kpix++)
65 {
66 jpix = kpix;
67 IX = 0;
68 IY = 0;
69 IP = 1 ;// ! bit position (in x and y)
70 while( jpix!=0 )
71 { // ! go through all the bits
72 ID=jpix%2;// ! bit value (in kpix), goes in ix
73 jpix = jpix/2;
74 IX = ID*IP+IX;
75
76 ID=jpix%2;// ! bit value (in kpix), goes in iy
77 jpix = jpix/2;
78 IY = ID*IP+IY;
79
80 IP = 2*IP;// ! next bit (in x and y)
81 }
82 pix2x_(kpix) = IX;// ! in 0,31
83 pix2y_(kpix) = IY;// ! in 0,31
84 }
85}
86
87void PIXELS_XY::mk_xy2pix()
88{
89 /*
90 =================================================
91 subroutine mk_xy2pix
92 =================================================
93 c sets the array giving the number of the pixel lying in (x,y)
94 c x and y are in {1,128}
95 c the pixel number is in {0,128**2-1}
96 c
97 c if i-1 = sum_p=0 b_p * 2^p
98 c then ix = sum_p=0 b_p * 4^p
99 c iy = 2*ix
100 c ix + iy in {0, 128**2 -1}
101 =================================================
102 */
103 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
104 // (16/12/98)
105
106 int K,IP,I,J,ID;
107 for(I = 1; I <= 128; I++)
108 {
109 J = I-1;// !pixel numbers
110 K = 0;//
111 IP = 1;//
112 truc : if( J==0 )
113 {
114 x2pix_(I-1) = K;
115 y2pix_(I-1) = 2*K;
116 }
117 else
118 {
119 ID = (int)fmod(J,2);
120 J = J/2;
121 K = IP*ID+K;
122 IP = IP*4;
123 goto truc;
124 }
125 }
126}
127
128
129
130int_4 nest2ring(int_4 nside, int_4 ipnest)
131{
132 /*
133 ====================================================
134 subroutine nest2ring(nside, ipnest, ipring)
135 ====================================================
136 c conversion from NESTED to RING pixel number
137 ====================================================
138 */
139 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
140 // (16/12/98)
141
142 const PIXELS_XY& PXY= PIXELS_XY::instance();
143
144 int npix, npface, face_num, ncap, n_before;
145 int ipf, ip_low, ip_trunc, ip_med, ip_hi;
146 int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
147 int ns_max=8192;
148 int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
149 int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
150
151 if( nside<1 || nside>ns_max ) {
152 cout << "nside out of range" << endl;
153 exit(0);
154 }
155 npix = 12 * nside* nside;
156 if( ipnest<0 || ipnest>npix-1 ) {
157 cout << "ipnest out of range" << endl;
158 exit(0);
159 }
160
161 ncap = 2* nside*( nside-1);// ! number of points in the North Polar cap
162 nl4 = 4* nside;
163
164 //c finds the face, and the number in the face
165 npface = nside* nside;
166 //cccccc ip = ipnest - 1 ! in {0,npix-1}
167
168 face_num = ipnest/npface;// ! face number in {0,11}
169 ipf =ipnest%npface;// ! pixel number in the face {0,npface-1}
170 //c finds the x,y on the face (starting from the lowest corner)
171 //c from the pixel number
172 ip_low=ipf%1024; // ! content of the last 10 bits
173 ip_trunc = ipf/1024; // ! truncation of the last 10 bits
174 ip_med=ip_trunc%1024; // ! content of the next 10 bits
175 ip_hi = ip_trunc/1024;// ! content of the high weight 10 bits
176
177 ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
178 iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
179
180 //c transforms this in (horizontal, vertical) coordinates
181 jrt = ix + iy;// ! 'vertical' in {0,2*(nside-1)}
182 jpt = ix - iy;// ! 'horizontal' in {-nside+1,nside-1}
183
184 //c computes the z coordinate on the sphere
185 // jr = jrll[face_num+1]*nside - jrt - 1;// ! ring number in {1,4*nside-1}
186 jr = jrll[face_num]*nside - jrt - 1;
187 nr = nside;// ! equatorial region (the most frequent)
188 n_before = ncap + nl4 * (jr - nside);
189 kshift=(jr - nside)%2;
190 if( jr<nside ) {//then ! north pole region
191 nr = jr;
192 n_before = 2 * nr * (nr - 1);
193 kshift = 0;
194 }
195 else if( jr>3*nside ) {//then ! south pole region
196 nr = nl4 - jr;
197 n_before = npix - 2 * (nr + 1) * nr;
198 kshift = 0;
199 }
200
201 //c computes the phi coordinate on the sphere, in [0,2Pi]
202 jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;// ! 'phi' number in the ring in {1,4*nr}
203
204 if( jp>nl4 ) jp = jp - nl4;
205 if( jp<1 ) jp = jp + nl4;
206
207 int aux=n_before + jp - 1;
208 return (n_before + jp - 1);// ! in {0, npix-1}
209}
210
211
212int_4 ring2nest(int_4 nside, int_4 ipring)
213{
214 /*
215 ==================================================
216 subroutine ring2nest(nside, ipring, ipnest)
217 ==================================================
218 c conversion from RING to NESTED pixel number
219 ==================================================
220 */
221 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
222 // (16/12/98)
223
224 const PIXELS_XY& PXY= PIXELS_XY::instance();
225
226 double fihip, hip;
227 int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
228 int kshift, face_num, nr;
229 int irn, ire, irm, irs, irt, ifm , ifp;
230 int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
231 int ns_max(8192);
232
233 // coordinate of the lowest corner of each face
234 int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
235 int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
236
237 if( nside<1 || nside>ns_max ) {
238 cout << "nside out of range" << endl;
239 exit(0);
240 }
241 npix = 12 * nside*nside;
242 if( ipring<0 || ipring>npix-1 ) {
243 cout << "ipring out of range" << endl;
244 exit(0);
245 }
246
247 nl2 = 2*nside;
248 nl4 = 4*nside;
249 npix = 12*nside*nside;// ! total number of points
250 ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
251 ipring1 = ipring + 1;
252
253 //c finds the ring number, the position of the ring and the face number
254 if( ipring1<=ncap ) {//then
255
256 hip = ipring1/2.;
257 fihip = floor ( hip );
258 irn = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
259 iphi = ipring1 - 2*irn*(irn - 1);
260
261 kshift = 0;
262 nr = irn ;// ! 1/4 of the number of points on the current ring
263 face_num = (iphi-1) / irn;// ! in {0,3}
264 }
265 else if( ipring1<=nl2*(5*nside+1) ) {//then
266
267 ip = ipring1 - ncap - 1;
268 irn = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
269 iphi = (int)fmod(ip,nl4) + 1;
270
271 kshift = (int)fmod(irn+nside,2);// ! 1 if irn+nside is odd, 0 otherwise
272 nr = nside;
273 ire = irn - nside + 1;// ! in {1, 2*nside +1}
274 irm = nl2 + 2 - ire;
275 ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
276 ifp = (iphi - irm/2 + nside -1) / nside;
277 if( ifp==ifm ) {//then ! faces 4 to 7
278 face_num = (int)fmod(ifp,4) + 4;
279 }
280 else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
281 face_num = ifp;
282 }
283 else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
284 face_num = ifp + 7;
285 }
286 }
287 else {
288
289 ip = npix - ipring1 + 1;
290 hip = ip/2.;
291 fihip = floor ( hip );
292 irs = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from South pole
293 iphi = 4*irs + 1 - (ip - 2*irs*(irs-1));
294
295 kshift = 0;
296 nr = irs;
297 irn = nl4 - irs;
298 face_num = (iphi-1) / irs + 8;// ! in {8,11}
299 }
300
301 //c finds the (x,y) on the face
302 irt = irn - jrll[face_num]*nside + 1;// ! in {-nside+1,0}
303 ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
304
305
306 if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
307
308 ix = (ipt - irt ) / 2;
309 iy = -(ipt + irt ) / 2;
310
311 ix_low = (int)fmod(ix,128);
312 ix_hi = ix/128;
313 iy_low = (int)fmod(iy,128);
314 iy_hi = iy/128;
315 ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
316
317 return (ipf + face_num* nside *nside);// ! in {0, 12*nside**2 - 1}
318}
319
320int_4 ang2pix_ring(int_4 nside, double theta, double phi)
321{
322 /*
323 ==================================================
324 c gives the pixel number ipix (RING)
325 c corresponding to angles theta and phi
326 c==================================================
327 */
328 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
329 // (16/12/98)
330
331 int nl2, nl4, ncap, npix, jp, jm, ipix1;
332 double z, za, tt, tp, tmp;
333 int ir, ip, kshift;
334
335 double piover2(Pi/2.);
336 double twopi(2.*Pi);
337 double z0(2./3.);
338 int ns_max(8192);
339
340 if( nside<1 || nside>ns_max ) {
341 cout << "nside out of range" << endl;
342 exit(0);
343 }
344
345 if( theta<0. || theta>Pi) {
346 cout << "theta out of range" << endl;
347 exit(0);
348 }
349
350 z = cos(theta);
351 za = fabs(z);
352 if( phi >= twopi) phi = phi - twopi;
353 if (phi < 0.) phi = phi + twopi;
354 tt = phi / piover2;// ! in [0,4)
355
356 nl2 = 2*nside;
357 nl4 = 4*nside;
358 ncap = nl2*(nside-1);// ! number of pixels in the north polar cap
359 npix = 12*nside*nside;
360
361 if( za <= z0 ) {
362
363 jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of ascending edge line
364 jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
365
366 ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
367 kshift = 0;
368 if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
369
370 ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
371 if( ip>nl4 ) ip = ip - nl4;
372
373 ipix1 = ncap + nl4*(ir-1) + ip ;
374 }
375 else {
376
377 tp = tt - floor(tt);// !MOD(tt,1.d0)
378 tmp = sqrt( 3.*(1. - za) );
379
380 jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
381 jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
382
383 ir = jp + jm + 1;// ! ring number counted from the closest pole
384 ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
385 if( ip>4*ir ) ip = ip - 4*ir;
386
387 ipix1 = 2*ir*(ir-1) + ip;
388 if( z<=0. ) {
389 ipix1 = npix - 2*ir*(ir+1) + ip;
390 }
391 }
392 return (ipix1 - 1);// ! in {0, npix-1}
393}
394
395int_4 ang2pix_nest(int_4 nside, double theta, double phi)
396{
397 /*
398 ==================================================
399 subroutine ang2pix_nest(nside, theta, phi, ipix)
400 ==================================================
401 c gives the pixel number ipix (NESTED)
402 c corresponding to angles theta and phi
403 c
404 c the computation is made to the highest resolution available (nside=8192)
405 c and then degraded to that required (by integer division)
406 c this doesn't cost more, and it makes sure
407 c that the treatement of round-off will be consistent
408 c for every resolution
409 ==================================================
410 */
411 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
412 // (16/12/98)
413
414 const PIXELS_XY& PXY= PIXELS_XY::instance();
415
416 double z, za, z0, tt, tp, tmp;
417 int face_num,jp,jm;
418 int ifp, ifm;
419 int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
420 double piover2(Pi/2.), twopi(2.*Pi);
421 int ns_max(8192);
422
423 if( nside<1 || nside>ns_max ) {
424 cout << "nside out of range" << endl;
425 exit(0);
426 }
427 if( theta<0 || theta>Pi ) {
428 cout << "theta out of range" << endl;
429 exit(0);
430 }
431 z = cos(theta);
432 za = fabs(z);
433 z0 = 2./3.;
434 if( phi>=twopi ) phi = phi - twopi;
435 if( phi<0. ) phi = phi + twopi;
436 tt = phi / piover2;// ! in [0,4[
437 if( za<=z0 ) { // then ! equatorial region
438
439 //(the index of edge lines increase when the longitude=phi goes up)
440 jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// ! ascending edge line index
441 jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
442
443 //c finds the face
444 ifp = jp / ns_max;// ! in {0,4}
445 ifm = jm / ns_max;
446 if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then ! faces 4 to 7
447 else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
448 else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
449
450 ix = (int)fmod(jm, ns_max);
451 iy = ns_max - (int)fmod(jp, ns_max) - 1;
452 }
453 else { //! polar region, za > 2/3
454
455 ntt = (int)floor(tt);
456 if( ntt>=4 ) ntt = 3;
457 tp = tt - ntt;
458 tmp = sqrt( 3.*(1. - za) );// ! in ]0,1]
459
460 //(the index of edge lines increase when distance from the closest pole goes up)
461 jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
462 jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
463 jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
464 jm = (int)min(ns_max-1, jm);
465
466 // finds the face and pixel's (x,y)
467 if( z>=0 ) {
468 face_num = ntt;// ! in {0,3}
469 ix = ns_max - jm - 1;
470 iy = ns_max - jp - 1;
471 }
472 else {
473 face_num = ntt + 8;// ! in {8,11}
474 ix = jp;
475 iy = jm;
476 }
477 }
478
479 ix_low = (int)fmod(ix,128);
480 ix_hi = ix/128;
481 iy_low = (int)fmod(iy,128);
482 iy_hi = iy/128;
483 ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
484 // ipf = ipf / pow(ns_max/nside,2.);// ! in {0, nside**2 - 1}
485 // return ( ipf + face_num*pow(nside,2));// ! in {0, 12*nside**2 - 1}
486 // $CHECK$ Reza 25/10/99 , pow remplace par *
487 ipf = ipf / ((ns_max/nside)*(ns_max/nside));
488 return (ipf + face_num*nside*nside);
489}
490
491void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) {
492 /*
493 ===================================================
494 c gives theta and phi corresponding to pixel ipix (RING)
495 c for a parameter nside
496 ===================================================
497 */
498 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
499 // (16/12/98)
500
501 int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
502 double fact1, fact2, fodd, hip, fihip;
503
504 int ns_max(8192);
505
506 if( nside<1 || nside>ns_max ) {
507 cout << "nside out of range" << endl;
508 exit(0);
509 }
510 npix = 12*nside*nside; // ! total number of points
511 if( ipix<0 || ipix>npix-1 ) {
512 cout << "ipix out of range" << endl;
513 exit(0);
514 }
515
516 ipix1 = ipix + 1; // in {1, npix}
517 nl2 = 2*nside;
518 nl4 = 4*nside;
519 ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
520 fact1 = 1.5*nside;
521 fact2 = 3.0*nside*nside;
522
523 if( ipix1 <= ncap ) { //! North Polar cap -------------
524
525 hip = ipix1/2.;
526 fihip = floor(hip);
527 iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
528 iphi = ipix1 - 2*iring*(iring - 1);
529
530 theta = acos( 1. - iring*iring / fact2 );
531 phi = ((double)iphi - 0.5) * Pi/(2.*iring);
532 // cout << theta << " " << phi << endl;
533 }
534 else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
535
536 ip = ipix1 - ncap - 1;
537 iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
538 iphi = ip%nl4 + 1;
539
540 fodd = 0.5 * (1 + (iring+nside)%2 );// ! 1 if iring+nside is odd, 1/2 otherwise
541 theta = acos( (nl2 - iring) / fact1 );
542 phi = ((double)iphi - fodd) * Pi /(2.*nside);
543 }
544 else {//! South Polar cap -----------------------------------
545
546 ip = npix - ipix1 + 1;
547 hip = ip/2.;
548 fihip = floor(hip);
549 iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from South pole
550 iphi = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
551
552 theta = acos( -1. + iring*iring / fact2 );
553 phi = ((double)iphi - 0.5) * Pi/(2.*iring);
554 // cout << theta << " " << phi << endl;
555 }
556}
557
558void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) {
559 /*
560 ==================================================
561 subroutine pix2ang_nest(nside, ipix, theta, phi)
562 ==================================================
563 c gives theta and phi corresponding to pixel ipix (NESTED)
564 c for a parameter nside
565 ==================================================
566 */
567 // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
568 // (16/12/98)
569
570 const PIXELS_XY& PXY= PIXELS_XY::instance();
571
572 int npix, npface, face_num;
573 int ipf, ip_low, ip_trunc, ip_med, ip_hi;
574 int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
575 double z, fn, fact1, fact2;
576 double piover2(Pi/2.);
577 int ns_max(8192);
578
579 // ! coordinate of the lowest corner of each face
580 int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
581 int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
582
583 if( nside<1 || nside>ns_max ) {
584 cout << "nside out of range" << endl;
585 exit(0);
586 }
587 npix = 12 * nside*nside;
588 if( ipix<0 || ipix>npix-1 ) {
589 cout << "ipix out of range" << endl;
590 exit(0);
591 }
592
593 fn = 1.*nside;
594 fact1 = 1./(3.*fn*fn);
595 fact2 = 2./(3.*fn);
596 nl4 = 4*nside;
597
598 //c finds the face, and the number in the face
599 npface = nside*nside;
600
601 face_num = ipix/npface;// ! face number in {0,11}
602 ipf = (int)fmod(ipix,npface);// ! pixel number in the face {0,npface-1}
603
604 //c finds the x,y on the face (starting from the lowest corner)
605 //c from the pixel number
606 ip_low = (int)fmod(ipf,1024);// ! content of the last 10 bits
607 ip_trunc = ipf/1024 ;// ! truncation of the last 10 bits
608 ip_med = (int)fmod(ip_trunc,1024);// ! content of the next 10 bits
609 ip_hi = ip_trunc/1024 ;//! content of the high weight 10 bits
610
611 ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
612 iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
613
614 //c transforms this in (horizontal, vertical) coordinates
615 jrt = ix + iy;// ! 'vertical' in {0,2*(nside-1)}
616 jpt = ix - iy;// ! 'horizontal' in {-nside+1,nside-1}
617
618 //c computes the z coordinate on the sphere
619 // jr = jrll[face_num+1]*nside - jrt - 1;// ! ring number in {1,4*nside-1}
620 jr = jrll[face_num]*nside - jrt - 1;
621 nr = nside;// ! equatorial region (the most frequent)
622 z = (2*nside-jr)*fact2;
623 kshift = (int)fmod(jr - nside, 2);
624 if( jr<nside ) { //then ! north pole region
625 nr = jr;
626 z = 1. - nr*nr*fact1;
627 kshift = 0;
628 }
629 else {
630 if( jr>3*nside ) {// then ! south pole region
631 nr = nl4 - jr;
632 z = - 1. + nr*nr*fact1;
633 kshift = 0;
634 }
635 }
636 theta = acos(z);
637
638 //c computes the phi coordinate on the sphere, in [0,2Pi]
639 // jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;// ! 'phi' number in the ring in {1,4*nr}
640 jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
641 if( jp>nl4 ) jp = jp - nl4;
642 if( jp<1 ) jp = jp + nl4;
643 phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
644}
645
646}
Note: See TracBrowser for help on using the repository browser.