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; -----------------------------------------------------------------------------
28pro test_tk2, nside, upix, random=random
29;+
30; test_tk2 [,nside ,upix, random= ]
31;
32;  tests the consistency of Healpix toolkit routines (pix2*, ang2*, vec2*, ang2vec, vec2ang)
33;
34;  by default, all pixels in 0, Npix-1 are tested (Npix=12*nside*nside)
35;  unless upix or random are defined
36;
37;  OPTIONAL INPUT
38;   nside (default = 32)
39;   upix , integer scalar or vector, list of hand-picked pixels
40;
41;  KEYWORD
42;   random : float scalar, some pixels (random*Npix +1) are picked randomly in [0, Npix-1]
43;
44;  EXAMPLES:
45;   test_tk2
46;   test_tk2, 1024
47;   test_tk2, 8192,  random=1.d-3
48;   test_tk2, 2L^29, random=1.d-12
49;
50;  HISTORY:
51;   based on test_tk, with a few extra tests on vectors
52;  2011-08: cosmetic editions
53;  2017-06: more cosmetic editions
54;-
55
56t0 = systime(1)
57if undefined(nside) then nside = 32
58;lnside = long(nside)
59
60trigger = (machar(/double)).eps * 8 ; 1.77e-15
61
62;npix = nside2npix(lnside,err=err_nside)
63npix = nside2npix(nside,err=err_nside)
64snpix = strtrim(string(npix,form='(i20)'),2)
65snpix1 = strtrim(string(npix-1,form='(i20)'),2)
66snside = strtrim(string(nside,form='(i9)'),2)
67
68if (err_nside ne 0) then begin
69    print,'Invalid Nside'
70    return
71endif
72
73if defined(random) then begin
74    nr = min([npix, long(random*npix)+1])
75    print,'Nside = '+snside
76    print,nr,' pixels are picked randomly in [0, '+snpix1+']'
77    if (nside gt 8192) then begin
78        pixel =long64( randomu(seed,nr, /double) * npix )
79    endif else begin
80        pixel =long( randomu(seed,nr, /double) * npix )
81    endelse
82    print,min(pixel),max(pixel)
83endif
84if defined(upix) then begin
85    if min(upix) ge 0 and max(upix) lt npix then begin
86        pixel = upix
87        print,'Nside = '+snside
88        print,'test selected pixels in ',min(upix), max(upix)
89    endif else begin
90        print,'invalid choice of pixels'
91        return
92    endelse
93endif
94if undefined(pixel) then begin
95    pixel = lindgen(npix)
96    print,'Nside = '+snside
97    print,'test all pixels in [0, '+snpix1+']'
98endif
99np = n_elements(pixel)
100
101error = 0
102
103;--------------------------
104pix2ang_ring, nside, pixel, theta, phi
105ang2pix_ring, nside, theta, phi, pixel2
106if total(abs(pixel2-pixel)) ne 0 then begin
107    print,'error pix <-> ang ring', nside
108    print,'out pixel range: ',minmax(pixel2)
109    error = error + 1
110endif
111
112pix2ang_nest, nside, pixel, theta, phi
113ang2pix_nest, nside, theta, phi, pixel2
114if total(abs(pixel2-pixel)) ne 0 then begin
115    print,'error pix <-> ang nest', nside
116    print,'out pixel range: ',minmax(pixel2)
117    error = error + 1
118endif
119;---------------------
120
121pix2vec_ring, nside, pixel, vec
122vec2pix_ring, nside, vec, pixel2
123if total(abs(pixel2-pixel)) ne 0 then begin
124    print,'error pix <-> vec ring', nside
125    error = error + 1
126    help,pixel, pixel2
127    bad = where( (pixel2-pixel) ne 0, nbad)
128    print,'#bad:',nbad
129    print,'in:',pixel[bad]
130    print,'out:',pixel2[bad]
131    print,'diff:',(pixel2-pixel)[bad]
132endif
133
134pix2vec_nest, nside, pixel, vec
135vec2pix_nest, nside, vec, pixel2
136if total(abs(pixel2-pixel)) ne 0 then begin
137    print,'error pix <-> vec nest', nside
138    print,'out pixel range: ',minmax(pixel2)
139    error = error + 1
140    bad = where( (pixel2-pixel) ne 0, nbad)
141    print,'#bad:',nbad
142;     print,pixel
143;     print,pixel-pixel2
144endif
145;---------------------
146
147ring2nest, nside, pixel, pix_n
148nest2ring, nside, pix_n, pixel2
149if total(abs(pixel2-pixel)) ne 0 then begin
150    print,'error ring <-> nest', nside
151    error = error + 1
152endif
153;---------------------
154
155
156pix2ang_ring, nside, pixel, theta, phi
157ang2vec,             theta, phi, vec
158vec2pix_nest, nside, vec, pixel1
159nest2ring,    nside, pixel1, pixel2
160if total(abs(pixel2-pixel)) ne 0 then begin
161    ;print,'error loop1', nside
162    print,'R -> A -> V -> N -> R', nside
163    error = error + 1
164endif
165;---------------------
166
167pix2vec_ring, nside, pixel, vec
168vec2ang,             vec, theta, phi
169ang2pix_nest, nside, theta, phi, pixel1
170nest2ring,    nside, pixel1, pixel2
171if total(abs(pixel2-pixel)) ne 0 then begin
172    ;print,'error loop2', nside
173    print,'R -> V -> A -> N -> R', nside
174    error = error + 1
175endif
176;---------------------
177
178pix2vec_ring, nside, pixel, vec, vv
179ring2nest, nside, pixel, pixnest
180pix2vec_nest, nside, pixnest, vec2, vv2
181t1 = total(abs(vec-vec2))/np
182m1 = max(abs(vec-vec2))
183if  (t1 gt trigger || m1 gt trigger) then begin
184    ;print,'error loop3 vector', nside
185    print,'R -> V  vs  R -> N -> V,  vector', nside
186    print,t1,m1,n_elements(vec)
187    print,total(vec,1)/np
188    error = error + 1
189endif
190t2 = total(abs(vv-vv2))/np
191m2 = max(abs(vv-vv2))
192if (t2 gt trigger || m2 gt trigger) then begin
193    ;print,'error loop3 vertex', nside
194    print,'R -> V  vs  R -> N -> V,  vertex', nside
195    print,t2,m2,n_elements(vv)
196    print,total(total(vv,1),2)/np
197    error = error + 1
198endif
199
200t1 = systime(1)
201if (error eq 0) then begin
202    print, ' -----------------'
203    print, '   TEST2 PASSED ('+snside+')'
204    print, ' -----------------'
205    print, 'Time [s]: ', t1-t0
206endif
207
208
209return
210end
211