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 data2azeq, data, pol_data, pix_type, pix_param, do_conv, do_rot, coord_in, coord_out, eul_mat, $ 29 color, Tmax, Tmin, color_bar, dx, planvec, vector_scale, $ 30 PXSIZE=pxsize, PYSIZE=pysize, ROT=rot_ang, LOG=log, HIST_EQUAL=hist_equal, $ 31 MAX=max_set, MIN=min_set, $ 32 RESO_ARCMIN=reso_arcmin, FITS = fits, $ 33 FLIP=flip, DATA_plot = data_plot, $ 34 POLARIZATION=polarization, SILENT=silent, PIXEL_LIST=pixel_list, ASINH=asinh, $ 35 TRUECOLORS=truecolors, DATA_TC=data_tc, MAP_OUT=map_out, HALF_SKY=half_sky 36 37;+ 38;============================================================================================== 39; DATA2AZEQ 40; 41; turns a Healpix or Quad-cube tessellation of the sphere 42; into a rectangular map in azimuthal equidistant coordinates 43; 44; DATA2AZEQ, data, pix_type, pix_param, do_conv, do_rot, coord_in, coord_out, eul_mat, 45; color, Tmax, Tmin, color_bar, dx, planvec, vector_scale, 46; pxsize=, pysize=, rot=, log=, hist_equal=, max=, min=, 47; reso_arcmin=, fits=, flip=, data_plot=, POLARIZATION=, SILENT=, 48; PIXEL_LIST=, ASINH=, TRUECOLORS=, DATA_TC=, MAP_OUT=, HALF_SKY= 49; 50; IN : 51; data, pix_type, pix_param, do_conv, do_rot, coord_in, coord_out, eul_mat 52; OUT : 53; color, Tmax, Tmin, color_bar, dx, planvec, vector_scale 54; KEYWORDS 55; Pxsize, Pysize, Rot, Log, Hist_equal, Max, Min, Reso_arcmin, 56; Fits, flip, data_plot, polarization, silent, pixel_list, asinh 57; 58; called by azeqview 59; 60; HISTORY 61; hacked from cartview 62;============================================================================================== 63;- 64 65do_true = keyword_set(truecolors) 66truetype = do_true ? truecolors : 0 67;help,data 68proj_small = 'azimequid' 69du_dv = 1. ; aspect ratio 70fudge = 1.00 ; 71if keyword_set(flip) then flipconv=1 else flipconv = -1 ; longitude increase leftward by default (astro convention) 72if undefined(polarization) then polarization=0 73do_polamplitude = (polarization[0] eq 1) 74do_poldirection = (polarization[0] eq 2) 75do_polvector = (polarization[0] eq 3) 76 77!P.BACKGROUND = 1 ; white background 78!P.COLOR = 0 ; black foreground 79 80mode_col = keyword_set(hist_equal) 81mode_col = mode_col + 2*keyword_set(log) + 4*keyword_value(asinh, default=0, min=0, max=2) 82 83obs_npix = N_ELEMENTS(data) 84npix_full = (pix_type eq 'Q') ? 6*(4L)^(pix_param-1) : nside2npix(pix_param) 85 86bad_data= !healpix.bad_value 87 88if (do_poldirection or do_polvector) then begin 89 ; compute new position of pixelisation North Pole in the plot coordinates 90 north_pole = [0.,0.,1.] 91 if (do_conv) then north_pole = SKYCONV(north_pole, inco= coord_in, outco=coord_out) 92 if (do_rot) then north_pole = north_pole # transpose(eul_mat) 93endif 94; ------------------------------------------------------------- 95; create the rectangular window 96; ------------------------------------------------------------- 97if defined(pxsize) then xsize = pxsize*1L else xsize = 500L 98if defined(pysize) then ysize = pysize*1L else ysize = xsize 99if defined(reso_arcmin) then resgrid = reso_arcmin/60. else resgrid = 1.5/60. 100dx = resgrid * !DtoR 101zsize = (do_true) ? 3 : 1 102N_uv = xsize*ysize 103indlist = (n_elements(pixel_list) eq n_elements(data[*,0])) 104dtype = size(data,/type) eq 5 ? 5 : 4 ; double (5) or float (4) by default 105 106if (~keyword_set(silent)) then begin 107 print,'Input map : ',3600.*6./sqrt(!dpi*npix_full),' arcmin / pixel ',form='(a,f8.3,a)' 108 print,'AzimEquid map :',resgrid*60.,' arcmin / pixel ',xsize,'*',ysize,form='(a,f8.3,a,i4,a,i4)' 109endif 110 111;;grid = FLTARR(xsize, ysize, zsize) 112;; grid = MAKE_ARRAY(/FLOAT,xsize,ysize, Value = bad_data) 113grid = MAKE_ARRAY(type=dtype, xsize, ysize, zsize) 114if do_polvector then planvec = MAKE_ARRAY(type=dtype,xsize,ysize, 2, Value = bad_data) 115; ------------------------------------------------------------- 116; makes the projection around the chosen contact point 117; ------------------------------------------------------------- 118; position on the planar grid 119x0 = +1. 120xll= 0 & xur = xsize-1 121yll= 0 & yur = ysize-1 122xc = 0.5*(xll+xur) 123yc = 0.5*(yll+yur) 124thetamax = keyword_set(half_sky) ? !PI*0.5 : !PI 125 126yband = LONG(5.e5 / FLOAT(xsize)) 127for ystart = 0, ysize - 1, yband do begin 128 yend = (ystart + yband - 1) < (ysize - 1) 129 nband = yend - ystart + 1 130 npb = xsize * nband 131 u = flipconv*(FINDGEN(xsize) - xc)# REPLICATE(dx,nband) ; minus sign = astro convention 132 v = REPLICATE(dx,xsize) # (FINDGEN(nband) + ystart - yc) 133 theta = sqrt(u*u + v*v) 134 off_mask = WHERE( theta gt thetamax, noff_mask) 135 if (noff_mask gt 0) then begin 136 if (undefined(plan_off)) then begin 137 plan_off = ystart*xsize+off_mask 138 endif else begin 139 plan_off = [plan_off, ystart*xsize+off_mask] 140 endelse 141 endif 142 cth = cos(theta) 143 sth = sin(theta) 144 x = reform(sth * u / (theta > 1.e-7), npb) 145 y = reform(sth * v / (theta > 1.e-7), npb) 146 z = reform(cth, npb) 147 ;vector = [[x],[y],[z]] ; normalised vector 148 vector = [[z],[x],[y]] ; normalised vector 149 ; -------------------------------- 150 ; deal with polarisation direction 151 ; -------------------------------- 152 if (do_poldirection or do_polvector) then begin 153 phi = 0. 154 if (do_rot or do_conv) then begin 155 vector = vector / (sqrt(total(vector^2, 2))#replicate(1,3)) ; normalize vector 156 ; compute rotation of local coordinates around each vector 157 tmp_sin = north_pole[1] * vector[*,0] - north_pole[0] * vector[*,1] 158 tmp_cos = north_pole[2] - vector[*,2] * (north_pole[0:2] ## vector) 159 if (flipconv lt 0) then tmp_cos = flipconv * tmp_cos 160 phi = ATAN(tmp_sin, tmp_cos) ; angle in radians 161 tmp_sin = 0. & tmp_cos = 0 162 endif 163 endif 164 ; --------- 165 ; rotation 166 ; --------- 167 if (do_rot) then vector = vector # eul_mat 168 if (do_conv) then vector = SKYCONV(vector, inco = coord_out, outco = coord_in) 169 ; we go from the final azeq map (system coord_out) to 170 ; the original one (system coord_in) 171 ; ----------------x--------------------------------------------- 172 ; converts the position on the sphere into pixel number 173 ; and project the corresponding data value on the map 174 ; ------------------------------------------------------------- 175 case pix_type of 176 'R' : VEC2PIX_RING, pix_param, vector, id_pix ; Healpix ring 177 'N' : VEC2PIX_NEST, pix_param, vector, id_pix ; Healpix nest 178 'Q' : id_pix = UV2PIX(vector, pix_param) ; QuadCube (COBE cgis software) 179 else : print,'error on pix_type' 180 endcase 181 if (do_true) then begin 182 for i=0,zsize-1 do grid[ystart*xsize+i*n_uv] = data_tc[id_pix,i] 183 endif else begin 184 if (do_poldirection) then begin 185 grid[ystart*xsize] = (data[id_pix] - phi + 4*!PI) MOD (2*!PI) ; in 0,2pi 186 endif else if (do_polvector) then begin 187 grid[ystart*xsize] = data[id_pix] 188 planvec[ystart*xsize] = pol_data[id_pix,0] 189 planvec[ystart*xsize+n_uv] = (pol_data[id_pix,1] - phi + 4*!PI) MOD (2*!PI) ; angle 190 endif else begin 191 ;grid[ystart*xsize] = data[id_pix] 192 grid[ystart*xsize] = sample_sparse_array(data,id_pix,in_pix=pixel_list,default=!healpix.bad_value) 193 endelse 194 endelse 195endfor 196u = 0 & v = 0 & x = 0 & vector = 0 197 198; ------------------------------------------------------------- 199; Test for unobserved pixels 200; ------------------------------------------------------------- 201data_plot = temporary(data) 202pol_data = 0 203find_min_max_valid, grid, mindata, maxdata, valid=Obs, bad_data=0.9*bad_data 204 205;----------------------------------- 206; export in FITS and as an array the original AzEq map before alteration 207;---------------------------------------------- 208 209; grid -> IDL array 210if arg_present(map_out) then map_out = proj2map_out(grid, offmap=plan_off, bad_data=bad_data) 211 212; grid -> FITS file 213if keyword_set(fits) then begin 214; message_patch,level=-1,'FITS not supported yet',/info 215 proj2fits, grid, fits, $ 216 projection = 'AZEQ', flip=flip, $ 217 rot = rot_ang, coord=coord_out, reso = resgrid, unit = sunits, min=mindata, max = maxdata 218endif 219 220; ------------------------------------------------------------- 221; set min and max and computes the color scaling 222; ------------------------------------------------------------- 223if (do_poldirection) then begin 224 min_set = 0. 225 max_set = 2*!pi 226endif 227 228if (truetype eq 2) then begin 229 ; truecolors=2 map each field to its color independently 230 color = bytarr(xsize,ysize,zsize) 231 for i=0,zsize-1 do begin 232 find_min_max_valid, grid[*,*,i], mindata, maxdata, valid=Obs, bad_data = 0.9 * bad_data 233 color[0,0,i] = COLOR_MAP(grid[*,*,i], mindata, maxdata, Obs, $ 234 color_bar = color_bar, mode=mode_col, silent=silent) 235 endfor 236endif else begin 237 ; same for truecolors=1 and false colors: 238 color = COLOR_MAP(grid, mindata, maxdata, Obs, $ 239 color_bar = color_bar, mode=mode_col, $ 240 minset = min_set, maxset = max_set, silent=silent) 241endelse 242 243if (defined(plan_off)) then begin 244 for i=0,zsize-1 do color[plan_off+i*n_uv] = !P.BACKGROUND ; white 245endif 246if (do_polvector) then begin ; rescale polarisation vector in each valid pixel 247 planvec[*,*,0] = vector_map(planvec[*,*,0], Obs, vector_scale = vector_scale) 248endif 249Obs = 0 250grid = 0 251Tmin = mindata & Tmax = maxdata 252 253return 254end 255 256