Source code for AegeanTools.angle_tools

#! /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 : string 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 : string 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 : string 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 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