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 ud_grade_cut4_old, pixel, signal, n_obs, serror, nside_in = nside_in, $
29          nside_out = nside_out, order_in = order_in, order_out = order_out, pessimistic = pessimistic
30;+
31;
32; UD_GRADE_CUT4_OLD,    $
33;     Pixel, Signal [, N_Obs, Serror],   $
34;     Nside_in=, [Nside_out=,] Order_in=, [Order_out=], [Pessimistic=]
35;
36;
37; up_gradation, de_gradation for the fields
38;  Pixel (pixel index),
39;  Signal (pixel temperature = average of several observations),
40;  N_obs (number of observations in pixel),
41;  Serror (error on temperature = rms of average)
42;
43;
44;  degradation :  K small pixels (i) -> 1 big pixel (j)
45;    N_Obs :    n_i  -> nf_j = sum_i (n_i)
46;    Signal:    t_i  -> tf_j = sum_i (n_i t_i)  / sum_i (n_i)
47;    Serror:
48;     e_i  -> ef_j = sqrt[ sum_i(n_i(t_i-tf_j)^2) + sum_i (n_i s_i)^2] / sum_i n_i
49;    Pixel : at least one exposed small pixel -> 1 exposed big pixel
50;
51;  upgradation : 1 big pixel(j) -> K small *identical* pixels (i)
52;    N_obs :    nf_j   -> n_i = nf_j / K (approximated to closest integer)
53;    Signal:    tf_j   -> t_i = tf_j
54;    Serror:    ef_j   -> e_i = ef_j * sqrt(K)
55;
56;
57;   EH, 2000-11-23
58;    June 2003,  EH, replaced STOPs by MESSAGEs
59;
60;    2009-03-25: fixed bug in prograding, allow large Nside
61;
62;-
63tstart = systime(1)
64routine = 'UD_GRADE_CUT4_OLD'
65if (undefined(nside_in) or undefined(order_in)) then begin
66    print,routine+': Pixel, Signal [, N_Obs, Serror], Nside_in=, [Nside_out=,] Order_in=, [Order_out=]'
67    print,'******** No processing ********'
68    return
69endif
70if undefined(nside_out) then nside_out = nside_in
71if undefined(order_out) then order_out = order_in
72npix_in = nside2npix(nside_in)
73npix_out = nside2npix(nside_out)
74if (npix_in le 0 or npix_out le 0) then begin
75    print,routine+': Invalid Nside_in or Nside_out :',nside_in,nside_out
76    message,'Abort'
77endif
78
79ord_in  = strmid(order_in,0,4)
80ord_out = strmid(order_out,0,4)
81
82if ((ord_in eq ord_out) and (nside_in eq nside_out)) then begin
83; ------------
84; nothing
85; ------------
86    return
87endif
88
89
90ok_obs = 0
91if defined(n_obs) then begin
92    if (max(n_obs) gt 0) then ok_obs = 1
93endif
94if (nside_out lt nside_in) then begin
95; ------------
96; degradation
97; ------------
98    ratio = (nside_in/nside_out)
99    ratio2 = ratio * ratio
100
101    if (ord_in eq 'RING') then begin
102        ring2nest, nside_in, pixel, pixel
103    endif
104
105    minhit = 1                                         ; at least one good sub-pixel
106    if (keyword_set(pessimistic)) then minhit = ratio2 ; all sub pixel should be good
107
108    hitin = bytarr(ratio2,npix_out)
109    hitin[pixel] = 1B
110    hitout = total(hitin,1)
111    if (ok_obs) then begin
112        obsin = lonarr(ratio2,npix_out)
113        obsin[pixel] = n_obs
114        obsout = total(obsin,1) ; total number of hits per large pixel
115    endif else begin
116        obsin = hitin
117        obsout = hitout         ; number of exposed subpixels per large pixel
118    endelse
119
120    mapin = fltarr(ratio2,npix_out)
121    mapin[pixel] = signal
122    mapout = total(obsin * mapin,1)
123
124    pixelout  = where(obsout gt 0 and hitout ge minhit, obs_npix)
125    signalout = mapout[pixelout]/obsout[pixelout] ; average temperature per large pixel
126
127    if defined(serror)  then begin
128        mapav = fltarr(ratio2,npix_out)
129        errin = fltarr(ratio2,npix_out)
130        ;;;mapav[pixel] = signalout[pixel/ratio2]
131        mapav[0:ratio2-1,pixelout] = replicate(1.0,ratio2)#signalout
132        errin[pixel] = serror
133        count = total(obsin gt 0, 1)>1.e-6 ; number of non-empty small pixels in big pixel
134        errout = sqrt(   total( obsin*errin*errin,1,/double)/count  + total(obsin*(mapin-mapav)^2,1,/double)     )
135        errin = 0
136        mapav = 0
137    endif
138
139    mapin = 0
140    obsin = 0
141
142    pixel = temporary(pixelout)
143    signal = temporary(signalout)
144    if defined(n_obs)  then n_obs = obsout[pixel]
145    if defined(serror) then serror = float(errout[pixel])/obsout[pixel]
146    nside = nside_out
147    npix = npix_out
148
149    if (ord_out eq 'RING') then begin
150        nest2ring, nside, pixel, pixr
151
152        ind = sort(pixr)
153        pixel = pixr[ind]
154        signal = signal[ind]
155        if defined(n_obs) then n_obs = n_obs[ind]
156        if defined(serror) then serror = serror[ind]
157        ind = 0
158    endif
159endif else begin
160; ------------
161; upgradation
162; ------------
163    ratio = (nside_out/nside_in)
164    ratio2 = ratio * ratio
165
166    ; go to NEST
167    if (ord_in eq 'RING') then begin
168        ring2nest, nside_in, pixel, pixel
169    endif
170
171    npix = n_elements(pixel)
172    npix_out = long64(npix) * ratio2
173
174    increm = lindgen(ratio2)
175    if ((nside_out le 8192) && (size(/type,pixel) ne 14)) then begin
176        pixel_out = lonarr(ratio2,npix) ; corrected 2009-03-25
177        lratio2   = ratio2
178    endif else begin
179        pixel_out = lon64arr(ratio2,npix)
180        lratio2   = long64(ratio2)
181    endelse
182    for i=0LL, npix-1LL do begin
183        pixel_out[*,i] = pixel[i] * lratio2 + increm
184    endfor
185    pixel  = reform(pixel_out, npix_out)
186    pixel_out = 0
187    increm = 0
188
189    subpix = replicate(1.,ratio2)
190    signal = signal ## (subpix)
191    signal = reform(signal, npix_out, /overwrite)
192
193    if defined(n_obs) then begin
194        n_obs  = n_obs  ## (subpix/ratio2)
195        n_obs  = reform(round(n_obs) , npix_out, /overwrite)
196    endif
197    if defined(serror) then begin
198        serror = serror ## (subpix*ratio)
199        serror = reform(serror, npix_out, /overwrite)
200    endif
201    npix = npix_out
202    nside = nside_out
203
204    ; go to output order
205    if (ord_out eq 'RING') then begin
206        nest2ring, nside, pixel, pixel
207    endif
208
209    ; sort
210    if (ord_in eq 'RING' or ord_out eq 'RING') then begin
211        ind = sort(pixel)
212        pixel = pixel[ind]
213        signal = signal[ind]
214        if defined(n_obs) then n_obs = n_obs[ind]
215        if defined(serror) then serror = serror[ind]
216    endif
217
218endelse
219print,'Time [s]:',systime(1)-tstart
220
221return
222end
223