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 cart2pix, xpos, ypos, id_pix, lon_deg, lat_deg, value
29;+
30; convert a uv position on a Cartesian map into a pixel number and
31; (long, lat)
32; only for scalar input
33;
34;-
35
36@viewcom ; define common
37do_poldirection = (polar_type eq 2)
38indlist = (n_elements(pixel_list) eq n_elements(data_plot))
39
40if (do_flip) then flipconv = +1 else flipconv = -1 ; longitude increase leftward by default (astro convention)
41
42if (do_poldirection) then begin
43    ; compute new position of pixelisation North Pole in the plot coordinates
44    north_pole = [0.,0.,1.]
45    if (do_conv) then north_pole = SKYCONV(north_pole, inco= coord_in, outco=coord_out)
46    if (do_rot) then north_pole = north_pole # transpose(eul_mat)
47endif
48
49id_pix  = -1
50lon_deg = -1000.
51lat_deg = -1000.
52
53if (xpos ge !x.crange(0) and xpos le !x.crange(1) and ypos ge !y.crange(0) and ypos le !y.crange(1)) then begin
54
55    vector = [cos(ypos)*cos(flipconv*xpos), cos(ypos)*sin(flipconv*xpos), sin(ypos)]
56    ; ---------
57    ; rotation
58    ; ---------
59    if (do_rot) then vector = vector # eul_mat
60    vec2ang, vector, lat_deg, lon_deg, /astro
61    if (do_conv) then vector = SKYCONV(vector, inco = coord_out, outco =  coord_in)
62          ; we go from the final Cartesian map (system coord_out) to
63          ; the original one (system coord_in)
64    ; -------------------------------------------------------------
65    ; converts the position on the sphere into pixel number
66    ; and project the corresponding data value on the map
67    ; -------------------------------------------------------------
68    case pix_type of
69        'R' : begin
70            vec2pix_ring, pix_param, vector, id_pix ; Healpix ring
71            pix2vec_ring, pix_param, id_pix, vectorp ; vector pointing to Healpix pixel center
72        end
73        'N' : begin
74            vec2pix_nest, pix_param, vector, id_pix ; Healpix nest
75            pix2vec_nest, pix_param, id_pix, vectorp ; vector pointing to Healpix pixel center
76        end
77        'Q' : begin
78            id_pix = uv2pix(vector, pix_param) ; QuadCube (COBE cgis software)
79            vectorp = pix2uv(id_pix, pix_param) ; vector pointing to QuadCube pixel center
80        end
81        else : print,'error on pix_type'
82    endcase
83    ; --------------------------------
84    ; deal with polarisation direction
85    ; --------------------------------
86    if (do_poldirection) then begin
87        phi = 0.
88        if (do_rot or do_conv) then begin
89            vector = vector / (sqrt(total(vector^2, 2))#replicate(1,3)) ; normalize vector
90            ; compute rotation of local coordinates around each vector
91            tmp_sin = north_pole[1] * vector[*,0] - north_pole[0] * vector[*,1]
92            tmp_cos = north_pole[2] - vector[*,2] * (north_pole[0:2] ## vector)
93            if (flipconv lt 0) then tmp_cos = flipconv * tmp_cos
94            phi = ATAN(tmp_sin, tmp_cos) ; angle in radians
95        endif
96        value = (data_plot[id_pix] - phi + 4*!PI) MOD (2*!PI) ; in 0,2pi
97    endif else begin
98        if (indlist) then begin
99            value = sample_sparse_array(data_plot, id_pix, in=pixel_list, default=!healpix.bad_value)
100        endif else begin
101            value = data_plot[id_pix]
102        endelse
103    endelse
104    if (n_elements(xpos) eq 1) then begin
105        id_pix = id_pix[0]
106        value = value[0]
107    endif
108endif
109
110
111return
112end
113
114