#!/usr/local/bin/python2.6 import pylab as pb import os import numpy as np import datetime import matplotlib ra = '17 45 40.045' dec = '-29 0 27.9' def utHoursToDateTime(ut): dt = [] for hr in ut: month = 3 # March day = 20 if (hr < 5.5): hr+=24 hr += 12 # We want the vernal equinox, not the autumnal while (hr >= 24): hr -= 24 day += 1 h = int(hr) m = int((hr-h)*60) s = int((hr-h-m/60.)*3600) mydate = datetime.datetime.strptime('%d-%d-2013 %d:%d:%d'%(day,month,h,m,s),'%d-%m-%Y %H:%M:%S') dt.append(mydate) return(dt) def mylabel(ut, el, name, color): pb.text(ut[list(el).index(np.max(el))-len(name)], np.max(el)+1, name, color=color, size=14) def process(fname): f = open(fname,'r') el = [] ut = [] lines = f.readlines() for line in lines: if (line.find('!') < 0): a,b,c,d,e,f,g,h,i = line.split() utc = float(b)+float(c)/60.+float(d)/3600. - hours ut.append(utc) el.append(float(h)) elif (line.find('longitude')>0): longitude = float(line.split('longitude =')[1]) print "longitude = ", longitude hours = longitude/15. return(np.array(el),np.array(ut)) obs = ['GLT','SMA','IRAM30m','PdBI','ALMA','LMT','CARMA','Haystack','SMT','SPT'] el = {} ut = {} time = {} pb.clf() desc = pb.subplot(111) for o in obs: cmd = 'verify -t -f %s.out -o %s -r %s -d %s'%(o,o,ra,dec) print "Running: %s" % (cmd) os.system(cmd) el[o],ut[o] = process('%s.out'%o) if (o == 'ALMA'): print ut[o] time[o] = pb.date2num(utHoursToDateTime(ut[o])) colors = ['b','g','k','k','r','m','y','c','b','c'] desc = pb.subplot(111) dashes = ['','','','-','','','','','','-'] pb.hold(True) for i in range(len(obs)): o = obs[i] if (o == 'SPT'): spt_el = abs(int(dec.split()[0])) pb.plot_date(pb.xlim(),[spt_el,spt_el],'%s-%s'%(colors[i],dashes[i])) else: pb.plot_date(time[o],el[o],'%s-%s'%(colors[i],dashes[i])) for i in range(len(obs)): o = obs[i] if (o == 'SPT'): pb.text(pb.xlim()[0]+1/24., spt_el+1, o, color=colors[i], size=14) else: try: mylabel(time[o], el[o], o, colors[i]) except: print "Failed on label for %s" % (o) desc.xaxis.set_major_locator(matplotlib.dates.HourLocator(byhour=range(0,24,2))) desc.xaxis.set_minor_locator(matplotlib.dates.HourLocator(byhour=range(0,24,1))) desc.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%H')) desc.xaxis.grid(True,which='major') desc.yaxis.grid(True,which='major') pb.xlabel('UT on Vernal Equinox') pb.ylabel('Elevation (deg)') pb.title('Sgr A*') png = 'sgra_eht.png' pb.savefig(png) print "Figure saved to %s" %png