#! /usr/bin/env python
"""
Tools for manipulating angles on the surface of a sphere
- distance
- bearing between two points
- translation along a path
- paths are either great circles or rhumb lines
also angle <-> string conversion tools for Aegean
"""
__author__ = "Paul Hancock"
import math
import numpy as np
[docs]def ra2dec(ra):
"""
Convert sexegessimal RA string into a float in degrees.
Parameters
----------
ra : str
A string separated representing the RA.
Expected format is `hh:mm[:ss.s]`
Colons can be replaced with any whit space character.
Returns
-------
ra : float
The RA in degrees.
"""
return dec2dec(ra)*15
[docs]def dec2dec(dec):
"""
Convert sexegessimal RA string into a float in degrees.
Parameters
----------
dec : str
A string separated representing the Dec.
Expected format is `[+- ]hh:mm[:ss.s]`
Colons can be replaced with any whit space character.
Returns
-------
dec : float
The Dec in degrees.
"""
d = dec.replace(':', ' ').split()
if len(d) == 2:
d.append('0.0')
if d[0].startswith('-') or float(d[0]) < 0:
return float(d[0]) - float(d[1]) / 60.0 - float(d[2]) / 3600.0
return float(d[0]) + float(d[1]) / 60.0 + float(d[2]) / 3600.0
[docs]def dec2dms(x):
"""
Convert decimal degrees into a sexagessimal string in degrees.
Parameters
----------
x : float
Angle in degrees
Returns
-------
dms : str
String of format [+-]DD:MM:SS.SS
or XX:XX:XX.XX if x is not finite.
"""
if not np.isfinite(x):
return 'XX:XX:XX.XX'
if x < 0:
sign = '-'
else:
sign = '+'
x = abs(x)
d = int(math.floor(x))
m = int(math.floor((x - d) * 60))
s = float(( (x - d) * 60 - m) * 60)
return '{0}{1:02d}:{2:02d}:{3:05.2f}'.format(sign, d, m, s)
[docs]def dec2hms(x):
"""
Convert decimal degrees into a sexagessimal string in hours.
Parameters
----------
x : float
Angle in degrees
Returns
-------
dms : string
String of format HH:MM:SS.SS
or XX:XX:XX.XX if x is not finite.
"""
if not np.isfinite(x):
return 'XX:XX:XX.XX'
# wrap negative RA's
if x < 0:
x += 360
x /= 15.0
h = int(x)
x = (x - h) * 60
m = int(x)
s = (x - m) * 60
return '{0:02d}:{1:02d}:{2:05.2f}'.format(h, m, s)
# The following functions are explained at http://www.movable-type.co.uk/scripts/latlong.html
# phi ~ lat ~ Dec
# lambda ~ lon ~ RA
[docs]def gcd(ra1, dec1, ra2, dec2):
"""
Calculate the great circle distance between to points using the haversine formula [1]_.
Parameters
----------
ra1, dec1, ra2, dec2 : float
The coordinates of the two points of interest.
Units are in degrees.
Returns
-------
dist : float
The distance between the two points in degrees.
Notes
-----
This duplicates the functionality of astropy but is faster as there is no creation of SkyCoords objects.
.. [1] `Haversine formula <https://en.wikipedia.org/wiki/Haversine_formula>`_
"""
# TODO: Vincenty formula see - https://en.wikipedia.org/wiki/Great-circle_distance
dlon = ra2 - ra1
dlat = dec2 - dec1
a = np.sin(np.radians(dlat) / 2) ** 2
a += np.cos(np.radians(dec1)) * np.cos(np.radians(dec2)) * np.sin(np.radians(dlon) / 2) ** 2
sep = np.degrees(2 * np.arcsin(np.minimum(1, np.sqrt(a))))
return sep
[docs]def bear(ra1, dec1, ra2, dec2):
"""
Calculate the bearing of point 2 from point 1 along a great circle.
The bearing is East of North and is in [0, 360), whereas position angle is also East of North but (-180,180]
Parameters
----------
ra1, dec1, ra2, dec2 : float
The sky coordinates (degrees) of the two points.
Returns
-------
bear : float
The bearing of point 2 from point 1 (degrees).
"""
rdec1 = np.radians(dec1)
rdec2 = np.radians(dec2)
rdlon = np.radians(ra2-ra1)
y = np.sin(rdlon) * np.cos(rdec2)
x = np.cos(rdec1) * np.sin(rdec2)
x -= np.sin(rdec1) * np.cos(rdec2) * np.cos(rdlon)
return np.degrees(np.arctan2(y, x))
[docs]def translate(ra, dec, r, theta):
"""
Translate a given point a distance r in the (initial) direction theta, along a great circle.
Parameters
----------
ra, dec : float
The initial point of interest (degrees).
r, theta : float
The distance and initial direction to translate (degrees).
Returns
-------
ra, dec : (float, float)
The translated position (degrees).
"""
factor = np.sin(np.radians(dec)) * np.cos(np.radians(r))
factor += np.cos(np.radians(dec)) * np.sin(np.radians(r)) * np.cos(np.radians(theta))
dec_out = np.degrees(np.arcsin(factor))
y = np.sin(np.radians(theta)) * np.sin(np.radians(r)) * np.cos(np.radians(dec))
x = np.cos(np.radians(r)) - np.sin(np.radians(dec)) * np.sin(np.radians(dec_out))
ra_out = ra + np.degrees(np.arctan2(y, x))
return ra_out, dec_out
[docs]def dist_rhumb(ra1, dec1, ra2, dec2):
"""
Calculate the Rhumb line distance between two points [1]_.
A Rhumb line between two points is one which follows a constant bearing.
Parameters
----------
ra1, dec1, ra2, dec2 : float
The position of the two points (degrees).
Returns
-------
dist : float
The distance between the two points along a line of constant bearing.
Notes
-----
.. [1] `Rhumb line <https://en.wikipedia.org/wiki/Rhumb_line>`_
"""
# verified against website to give correct results
phi1 = np.radians(dec1)
phi2 = np.radians(dec2)
dphi = phi2 - phi1
lambda1 = np.radians(ra1)
lambda2 = np.radians(ra2)
dpsi = np.log(np.tan(np.pi / 4 + phi2 / 2) / np.tan(np.pi / 4 + phi1 / 2))
if dpsi < 1e-12:
q = np.cos(phi1)
else:
q = dpsi / dphi
dlambda = lambda2 - lambda1
if dlambda > np.pi:
dlambda -= 2 * np.pi
dist = np.hypot(dphi, q * dlambda)
return np.degrees(dist)
[docs]def bear_rhumb(ra1, dec1, ra2, dec2):
"""
Calculate the bearing of point 2 from point 1 along a Rhumb line.
The bearing is East of North and is in [0, 360), whereas position angle is also East of North but (-180,180]
Parameters
----------
ra1, dec1, ra2, dec2 : float
The sky coordinates (degrees) of the two points.
Returns
-------
dist : float
The bearing of point 2 from point 1 along a Rhumb line (degrees).
"""
# verified against website to give correct results
phi1 = np.radians(dec1)
phi2 = np.radians(dec2)
lambda1 = np.radians(ra1)
lambda2 = np.radians(ra2)
dlambda = lambda2 - lambda1
dpsi = np.log(np.tan(np.pi / 4 + phi2 / 2) / np.tan(np.pi / 4 + phi1 / 2))
theta = np.arctan2(dlambda, dpsi)
return np.degrees(theta)
[docs]def translate_rhumb(ra, dec, r, theta):
"""
Translate a given point a distance r in the (initial) direction theta, along a Rhumb line.
Parameters
----------
ra, dec : float
The initial point of interest (degrees).
r, theta : float
The distance and initial direction to translate (degrees).
Returns
-------
ra, dec : float
The translated position (degrees).
"""
# verified against website to give correct results
# with the help of http://williams.best.vwh.net/avform.htm#Rhumb
delta = np.radians(r)
phi1 = np.radians(dec)
phi2 = phi1 + delta * np.cos(np.radians(theta))
dphi = phi2 - phi1
if abs(dphi) < 1e-9:
q = np.cos(phi1)
else:
dpsi = np.log(np.tan(np.pi / 4 + phi2 / 2) / np.tan(np.pi / 4 + phi1 / 2))
q = dphi / dpsi
lambda1 = np.radians(ra)
dlambda = delta * np.sin(np.radians(theta)) / q
lambda2 = lambda1 + dlambda
ra_out = np.degrees(lambda2)
dec_out = np.degrees(phi2)
return ra_out, dec_out