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 orth2pix, xpos, ypos, id_pix, lon_deg, lat_deg, value
29;+
30; convert a uv position on a Orthographic 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
53
54if (do_fullsky) then begin
55    c0 =  1
56    c1 = -1
57    valid = ((xpos-c0)^2 + ypos^2 le 1.) or ((xpos-c1)^2 + ypos^2 le 1.)
58endif else begin
59    valid = (xpos)^2 + ypos^2 le 1.
60endelse
61
62if (valid) then begin
63
64    nested = keyword_set(nest)
65    v1 = ypos
66    u1 = flipconv * xpos
67    if (do_fullsky) then begin
68        sign = (u1 ge 0.)*2 - 1 ; =1 for u1>0 , =-1 otherwise
69        ys = abs(u1)-c0
70        xs = sqrt(1.-ys*ys-v1*v1)
71        vector = [[sign*xs],[ys],[v1]] ; normalized vector
72    endif else begin
73        xs = sqrt(1.-u1*u1-v1*v1)
74        vector = [[xs],[u1],[v1]] ; normalized vector
75    endelse
76        ; ---------
77        ; rotation
78        ; ---------
79    if (do_rot) then vector = vector # eul_mat
80    vec2ang, vector, lat_deg, lon_deg, /astro
81    if (do_conv) then vector = SKYCONV(vector, inco = coord_out, outco =  coord_in)
82                                ; we go from the final Orthographic map (system coord_out) to
83                                ; the original one (system coord_in)
84        ; -------------------------------------------------------------
85        ; converts the position on the sphere into pixel number
86        ; and project the corresponding data value on the map
87        ; -------------------------------------------------------------
88    case pix_type of
89        'R' : begin
90            vec2pix_ring, pix_param, vector, id_pix ; Healpix ring
91            pix2vec_ring, pix_param, id_pix, vectorp ; vector pointing to Healpix pixel center
92        end
93        'N' : begin
94            vec2pix_nest, pix_param, vector, id_pix ; Healpix nest
95            pix2vec_nest, pix_param, id_pix, vectorp ; vector pointing to Healpix pixel center
96        end
97        'Q' : begin
98            id_pix = uv2pix(vector, pix_param) ; QuadCube (COBE cgis software)
99            vectorp = pix2uv(id_pix, pix_param) ; vector pointing to QuadCube pixel center
100        end
101        else : print,'error on pix_type'
102    endcase
103        ; --------------------------------
104        ; deal with polarisation direction
105        ; --------------------------------
106    if (do_poldirection) then begin
107        phi = 0.
108        if (do_rot or do_conv) then begin
109                                ; rotate pixel center to final coordinate system
110            if (do_conv) then vectorp = SKYCONV(vectorp, inco= coord_in, outco=coord_out)
111            if (do_rot) then vectorp = vectorp # transpose(eul_mat)
112                                ; compute rotation of local coordinates around each vector
113            tmp_sin = north_pole[1] * vectorp[*,0] - north_pole[0] * vectorp[*,1]
114            tmp_cos = north_pole[2] - vectorp[*,2] * (north_pole[0:2] ## vectorp)
115            if (flipconv lt 0) then tmp_cos = flipconv * tmp_cos
116            phi = ATAN(tmp_sin, tmp_cos) ; angle in radians
117        endif
118        value = (data_plot[id_pix] - phi + 4*!PI) MOD (2*!PI) ; in 0,2pi
119    endif else begin
120        if (indlist) then begin
121            value = sample_sparse_array(data_plot, id_pix, in=pixel_list, default=!healpix.bad_value)
122        endif else begin
123            value = data_plot[id_pix]
124        endelse
125    endelse
126    if (n_elements(xpos) eq 1) then begin
127        id_pix = id_pix[0]
128        value = value[0]
129    endif
130endif
131
132return
133end
134