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