1# This file is part of the Astrometry.net suite.
2# Licensed under a 3-clause BSD style license - see LICENSE
3
4from math import *  # Needed if we're going to exclude pylab :)
5
6PIl = 3.1415926535897932384626433832795029
7def rad2deg(r):    return 180.0*r/PIl
8def deg2rad(d):    return d*PIl/180.0
9def rad2arcmin(r): return 10800.0*r/PIl
10def arcmin2rad(a): return a*PIl/10800.0
11def rad2arcsec(r): return 648000.0*r/PIl
12def arcsec2rad(a): return a*PIl/648000.0
13def radec2x(r,d):  return cos(d)*cos(r) # r,d in radians
14def radec2y(r,d):  return cos(d)*sin(r) # r,d in radians
15def radec2z(r,d):  return sin(d)        # r,d in radians
16def z2dec(z):      return arcsin(z)     # result in radians
17def xy2ra(x,y):
18    "Convert x,y to ra in radians"
19    r = arctan2(y,x)
20    r += 2*PIl*(r>=0.00)
21    return r
22
23# wrapper to give us a tuple.
24def radec2xyz(r,d): return (radec2x(r,d), radec2y(r,d), radec2z(r,d))
25
26# Takes min/max RA & DEC (which a number of the catalogs provide) and
27# compute the xyz coordinates of the box on the surface of the sphere.
28def radecbox2xyz(min_r, max_r, min_d, max_d):
29    return [radec2xyz(min_r, min_d), radec2xyz(min_r, max_d), \
30        radec2xyz(max_r, max_d), radec2xyz(max_r, min_d)]
31
32
33
34#define radscale2xyzscale(r) (sqrt(2.0-2.0*cos(r/2.0)))
35def star_coords(s,r):
36    # eta is a vector perpendicular to r
37    etax = -r[1]
38    etay = +r[0]
39    etaz = 0.0
40    eta_norm = sqrt(etax * etax + etay * etay)
41    etax /= eta_norm
42    etay /= eta_norm
43
44    # xi =  r cross eta
45    xix = -r[2] * etay
46    xiy =  r[2] * etax
47    xiz =  r[0] * etay - r[1] * etax
48    sdotr = dot(s,r)
49
50    return (
51            s[0] * xix / sdotr +
52            s[1] * xiy / sdotr +
53            s[2] * xiz / sdotr,
54            s[0] * etax / sdotr +
55            s[1] * etay / sdotr
56           )
57
58def radectohealpix(ra, dec):
59    x = radec2x(ra, dec)
60    y = radec2y(ra, dec)
61    z = radec2z(ra, dec)
62    return xyztohealpix(x, y, z)
63
64def xyztohealpix(x, y, z) :
65    "Return the healpix catalog number associated with point (x,y,z)"
66
67    # the polar pixel center is at (z,phi/pi) = (2/3, 1/4)
68    # the left pixel center is at (z,phi/pi) = (0, 0)
69    # the right pixel center is at (z,phi/pi) = (0, 1/2)
70    twothirds = 2.0 / 3.0
71
72    phi = arctan2(y, x)
73    if phi < 0.0:
74        phi += 2.0 * pi
75
76    phioverpi = phi / pi
77
78    if z >= twothirds or z <= -twothirds:
79        # definitely in the polar regions.
80        # the north polar healpixes are 0,1,2,3
81        # the south polar healpixes are 8,9,10,11
82        if z >= twothirds:
83            offset = 0
84        else:
85            offset = 8
86
87        pix = int(phioverpi * 2.0)
88
89        return offset + pix
90
91    # could be polar or equatorial.
92    offset = int(phioverpi * 2.0)
93    phimod = phioverpi - offset * 0.5
94
95    z1 =  twothirds - (8.0 / 3.0) * phimod
96    z2 = -twothirds + (8.0 / 3.0) * phimod
97
98    if z >= z1 and z >= z2:
99        # north polar
100        return offset
101    if z <= z1 and z <= z2:
102        # south polar
103        return offset + 8
104    if phimod < 0.25:
105        # left equatorial
106        return offset + 4
107    # right equatorial
108    return ((offset+1)%4) + 4
109
110
111# Brought over from starutil.c & plotcat.c
112def project_equal_area(point):
113    x, y, z = point
114    xp = x * sqrt(1. / (1. + z))
115    yp  = y * sqrt(1. / (1. + z))
116    xp = 0.5 * (1.0 + xp)
117    yp = 0.5 * (1.0 + yp)
118    return (xp, yp)
119
120def project_hammer_aitoff_x(point):
121    x, y, z = point
122    theta = atan(float(x) / (z + 0.000001))
123    r = sqrt(x * x + z * z)
124    if z < 0:
125        if x < 0:
126            theta = theta - pi
127        else:
128            theta = pi + theta
129    theta /= 2.0
130    zp = r * cos(theta)
131    xp = r * sin(theta)
132    assert zp >= -0.01
133    return project_equal_area((xp, y, zp))
134
135
136def getxy(projectedpoint, N):
137    px, py = projectedpoint
138    px = 0.5 + (px - 0.5) * 0.99
139    py = 0.5 + (py - 0.5) * 0.99
140    X = round(px * 2*N)
141    Y = round(py * N)
142    return (X, Y)
143
144