diff options
Diffstat (limited to 'python/radec_angle.py')
| -rwxr-xr-x | python/radec_angle.py | 248 |
1 files changed, 0 insertions, 248 deletions
diff --git a/python/radec_angle.py b/python/radec_angle.py deleted file mode 100755 index ec01807..0000000 --- a/python/radec_angle.py +++ /dev/null @@ -1,248 +0,0 @@ -#!/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() |