Add radec2deg.py and radec_angle.py

1 files changed, 248 insertions, 0 deletions

diff --git a/python/radec_angle.py b/python/radec_angle.py
new file mode 100755
index 0000000..ec01807
--- /dev/null
+++ b/python/radec_angle.py
@@ -0,0 +1,248 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# Created: 2015-04-17
+# Updated: 2016-06-30
+#
+
+"""
+Calculate the angles between the given center point to other points
+on the sphere.
+
+The following input formats are supported:
+    longitude  latitude     => FMT
+    ------------------------------------
+    ??h??m??s  ??d??m??s    => "radec"
+    ?? ?? ??   ?? ?? ??     => "ra3dec3"
+    deg        deg          => "degdeg"
+"""
+
+import os
+import sys
+import re
+import getopt
+import math
+
+
+USAGE = """Usage:
+    %(prog)s [ -b -h ] -r RA -d DEC -i infile -f FMT -u UNIT
+
+Required arguments:
+    -r, --ra
+        RA (??h??m??s) of the center
+    -d, --dec
+        DEC (??d??m??s) of the center
+    -i, --infile
+        input file containing the coordinates data
+    -f, --format
+        value: radec | ra3dec3 | degdeg
+        coordinates format of the input data file
+    -u, --unit
+        value: deg | arcmin | arcsec
+        unit of the output data
+
+Optional arguments:
+    -b, --brief
+        brief mode: only output results
+    -h, --help
+""" % {'prog': os.path.basename(sys.argv[0])}
+
+
+def usage():
+    print(USAGE, file=sys.stderr)
+
+
+def ra2deg(h, m, s):
+    """
+    Convert RA (hour, minute, second) to degree.
+    """
+    return h * 15.0 + m * 15.0/60.0 + s * 15.0/3600.0
+
+
+def dec2deg(d, m, s):
+    """
+    Convert DEC (deg, arcmin, arcsec) to degree.
+    """
+    if (d >= 0):
+        sign = 1.0
+    else:
+        sign = -1.0
+    return sign * (math.fabs(d) + m/60.0 + s/3600.0)
+
+
+def s_ra2deg(hms):
+    """
+    Convert RA string ("??h??m??s") to degree.
+    """
+    h, m, s = map(float, re.sub('[hms]', ' ', hms).split())
+    return ra2deg(h, m, s)
+
+
+def s_dec2deg(dms):
+    """
+    Convert DEC string ("??d??m??s") to degree.
+    """
+    d, m, s = map(float, re.sub('[dms]', ' ', dms).split())
+    return dec2deg(d, m, s)
+
+
+def deg2rad(deg):
+    """
+    Convert unit from deg to rad.
+    """
+    return deg * math.pi / 180.0
+
+
+def rad2deg(rad):
+    """
+    Convert unit from rad to deg.
+    """
+    return rad * 180.0 / math.pi
+
+
+def central_angle(p1, p2, unit="deg"):
+    """
+    Calculate the central angle between the two points on the sphere.
+
+    Input parameters:
+        p1, p2: (longitude, latitude), coorindates of the two points
+
+    Algorithm:
+        (radial, azimuthal, polar): (r, theta, phi)
+        central_angle: alpha
+        longitude: lambda = theta
+        latitude: delta = 90 - phi
+        colatitude: phi
+
+        Unit vector:
+        \hat{r}_1 = (cos(theta1) sin(phi1), sin(theta1) sin(phi1), cos(phi1))
+            = (cos(lambda1) cos(delta1), sin(lambda1) cos(delta1), sin(delta1))
+        \hat{r}_2 = (cos(theta2) sin(phi2), sin(theta2) sin(phi2), cos(phi2))
+            = (cos(lambda2) cos(delta2), sin(lambda2) cos(delta2), sin(delta2))
+
+        Therefore the angle (alpha) between \hat{r}_1 and \hat{r}_2:
+        cos(alpha) = \hat{r}_1 \cdot \hat{r}_2
+            = cos(delta1) cos(delta2) cos(lambda1-lambda2)
+              + sin(delta1) sin(delta2)
+
+    References:
+    [1] Spherical Coordinates - Wolfram MathWorld
+        http://mathworld.wolfram.com/SphericalCoordinates.html
+        Equation (19)
+    [2] Great Circle - Wolfram MathWorld
+        http://mathworld.wolfram.com/GreatCircle.html
+        Equation (1), (2), (4)
+    """
+    lbd1, delta1 = map(deg2rad, p1)
+    lbd2, delta2 = map(deg2rad, p2)
+    dotvalue = (math.cos(delta1) * math.cos(delta2) * math.cos(lbd1-lbd2) +
+                math.sin(delta1) * math.sin(delta2))
+    alpha = math.acos(dotvalue)
+    if unit == "rad":
+        return alpha
+    elif unit == "arcmin":
+        return rad2deg(alpha) * 60.0
+    elif unit == "arcsec":
+        return rad2deg(alpha) * 60.0*60.0
+    else:
+        # default: degree
+        return rad2deg(alpha)
+
+
+def main():
+    # Mandatory arguments
+    center_ra = None
+    center_dec = None
+    infile = None
+
+    # Default parameters
+    unit = "arcmin"
+    fmt = "radec"  # default format: "??h??m??s  ??d??m??s"
+
+    # Process command line arguments
+    try:
+        opts, args = getopt.getopt(sys.argv[1:], "bd:f:hi:r:u:",
+                                   ["brief", "dec=", "format=", "help",
+                                    "infile=", "ra=", "unit="])
+    except getopt.GetoptError as err:
+        print(err)
+        usage()
+        sys.exit(2)
+    brief = False  # brief mode
+    for opt, arg in opts:
+        if opt in ("-h", "--help"):
+            usage()
+            sys.exit(1)
+        elif opt in ("-b", "--brief"):
+            brief = True
+        elif opt in ("-d", "--dec"):
+            center_dec = arg
+        elif opt in ("-r", "--ra"):
+            center_ra = arg
+        elif opt in ("-i", "--infile"):
+            infile = arg
+        elif opt in ("-f", "--format"):
+            fmt = arg
+        elif opt in ("-u", "--unit"):
+            unit = arg
+        else:
+            assert False, "unhandled option"
+
+    # Check mandatory arguments
+    if not center_ra:
+        print("Error: --ra argument required!", file=sys.stderr)
+    if not center_dec:
+        print("Error: --dec argument required!", file=sys.stderr)
+    if not infile:
+        print("Error: --infile argument required!", file=sys.stderr)
+
+    if fmt == "radec":
+        center_ra_deg = s_ra2deg(center_ra)
+        center_dec_deg = s_dec2deg(center_dec)
+    elif fmt == "ra3dec3":
+        ra_h, ra_m, ra_s = map(float, center_ra.split())
+        dec_d, dec_m, dec_s = map(float, center_dec.split())
+        center_ra_deg = ra2deg(ra_h, ra_m, ra_s)
+        center_dec_deg = dec2deg(dec_d, dec_m, dec_s)
+    elif fmt == "degdeg":
+        center_ra_deg = float(center_ra)
+        center_dec_deg = float(center_dec)
+    else:
+        print("Error: unknown format type: %s" % fmt, file=sys.stderr)
+        sys.exit(2)
+
+    if not brief:
+        print("# Central_angle (unit: %s)" % unit)
+
+    datafile = open(infile, "r")
+    for line in datafile:
+        if re.match(r"^\s*#", line):
+            # skip comments
+            continue
+        elif re.match(r"^\s*$", line):
+            # skip blank line
+            continue
+        # coordinate format
+        if fmt == "radec":
+            ra, dec = line.split()
+            ra_deg = s_ra2deg(ra)
+            dec_deg = s_dec2deg(dec)
+        elif fmt == "ra3dec3":
+            ra_h, ra_m, ra_s, dec_d, dec_m, dec_s = map(float, line.split())
+            ra_deg = ra2deg(ra_h, ra_m, ra_s)
+            dec_deg = dec2deg(dec_d, dec_m, dec_s)
+        elif fmt == "degdeg":
+            ra_deg, dec_deg = map(float, line.split())
+        else:
+            print("Error: unknown format type: %s" % fmt, file=sys.stderr)
+            sys.exit(2)
+        # calculate angle
+        angle = central_angle((center_ra_deg, center_dec_deg),
+                              (ra_deg, dec_deg), unit=unit)
+        print("%.10f" % angle)
+    datafile.close()
+
+
+if __name__ == "__main__":
+    main()