diff options
Diffstat (limited to 'python/radec_angle.py')
| -rwxr-xr-x | python/radec_angle.py | 248 |
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() |