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