#!/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()