#!/usr/local/bin/python2.6 # This script takes the output of the C program azelrate.c and # finds the points with azimuth velocity less than the specified # threshold, then sorts the result so that it can be plotted. import os, math import pylab as pb from matplotlib.ticker import MultipleLocator threshold = 1/3600. # degree/minute infile = 'azelrate.out.sma' def angularSeparation(ra0,dec0,ra1,dec1, returnComponents=False): """ Usage: au.angularSeparation(ra0,dec0,ra1,dec1) Computes the great circle angle between two celestial coordinates. using the Vincenty formula (from wikipedia) which is correct for all angles, as long as you use atan2() to handle a zero denominator. See http://en.wikipedia.org/wiki/Great_circle_distance ra,dec must be given in degrees, as is the output. It also works for the az,el coordinate system. Comopnent separations are field_0 minus field_1. See also angularSeparationRadians() -- Todd Hunter """ ra0 *= math.pi/180. dec0 *= math.pi/180. ra1 *= math.pi/180. dec1 *= math.pi/180. deltaLong = ra0-ra1 argument1 = (((math.cos(dec1)*math.sin(deltaLong))**2) + ((math.cos(dec0)*math.sin(dec1)-math.sin(dec0)*math.cos(dec1)*math.cos(deltaLong))**2))**0.5 argument2 = math.sin(dec0)*math.sin(dec1) + math.cos(dec0)*math.cos(dec1)*math.cos(deltaLong) angle = math.atan2(argument1, argument2) / (math.pi/180.) if (returnComponents): radegrees = (ra0-ra1)*180/math.pi decdegrees = (dec0-dec1)*180/math.pi retval = angle,radegrees,decdegrees else: retval = angle return(retval) lines = open(infile) outfile1 = infile + '.temp' outfile2 = infile + '.sort' out = open(outfile1,'w') for line in lines: if (line[0] == '#'): tokens = line.split() latitude = float(tokens[6]) continue a,b,c,d = line.split() az = float(a) el = float(b) angleToPole = angularSeparation(az,el,0,latitude) azrate = float(c) elrate = float(d) if (abs(azrate) < threshold): if (az > 180): az -=360 declination = 90 - angleToPole # print "%f %f %f" % (az,el,declination) out.write("%f %f %f\n" % (az,el,declination)) out.write("0.00 %f 90.00\n" % (latitude)) out.close() os.system('sort -n -k 1 %s > %s' % (outfile1,outfile2)) os.system('rm %s' % (outfile1)) print "result left in %s" % (outfile2) pb.clf() adesc = pb.subplot(111) lines = open(outfile2) az = [] el = [] declination = [] for line in lines: a,b,c = line.split() az.append(float(a)) el.append(float(b)) declination.append(float(c)) pb.plot(az,el,'b-',az,declination,'r-') pb.xlabel('Azimuth (deg)') pb.ylabel('Vertical angle (deg)') adesc.xaxis.grid(True,which='major') adesc.yaxis.grid(True,which='major') majorLocator = MultipleLocator(10) adesc.xaxis.set_major_locator(majorLocator) majorLocator = MultipleLocator(5) adesc.yaxis.set_major_locator(majorLocator) y0,y1 = pb.ylim() pb.ylim([y0,y1+1]) pb.text(0.1,0.95,'Elevation',color='b',transform=adesc.transAxes) pb.text(0.1,0.90,'Declination',color='r',transform=adesc.transAxes) pb.title('Directions with zero sidereal azimuth velocity at latitude %+f' % latitude) pb.savefig(outfile2+'.png', dpi=200) print "plot left in %s" % (outfile2+'.png') pb.draw()