#!/usr/local/bin/python2.6 import pylab as pb import os import numpy as np import datetime import matplotlib def utHoursToDateTime(ut): dt = [] for hr in ut: day = 20 hr += 12 if (hr >= 24): hr -= 24 day = 21 h = int(hr) m = int((hr-h)*60) s = int((hr-h-m/60.)*3600) mydate = datetime.datetime.strptime('%d-03-2013 %d:%d:%d'%(day,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 # utc += 12 # we want the vernal equinox # utc += 20 # we want Nov 20 # if (utc<0): # utc+=24 # if (utc>=24): # utc-=24 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)) ra = '12 30 49.4' dec = '12 23 28' obs = ['GLT','SMA','IRAM30m','PDBI','ALMA','LMT','CARMA','Haystack','SMT'] 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) time[o] = pb.date2num(utHoursToDateTime(ut[o])) colors = ['b','g','k','k','r','m','y','c','b'] desc = pb.subplot(111) dashes = ['','','','-','','','','','-'] pb.hold(True) for i in range(len(obs)): o = obs[i] pb.plot_date(time[o],el[o],'%s-%s'%(colors[i],dashes[i])) mylabel(time[o], el[o], o, colors[i]) 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('M87') png = 'm87_eht.png' pb.savefig(png) print "Figure saved to %s" %png