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