Skymaps now include start at injection location

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@65 c56465c9-8126-4a4f-9d7d-ac845eff4865

postprocess/skymap/skygram.py

+
 import numpy as np
 import healpy as hp
 import matplotlib
@@ -5,10 +6,12 @@ matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import matplotlib.mlab as mlab
 import find_trig as ft
+from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
+from pylal.xlal.date import XLALGreenwichSiderealTime
+
 
 plt.close('all')
-# NSIDE = 64
-NSIDE = 16
+NSIDE = 64
 print "The number of pixels is {0}".format(hp.nside2npix(NSIDE))
 
 # -------------
@@ -27,10 +30,27 @@ for theta, ra in zip(thetalist, ralist):
 hp.mollview(m,title="Sweet skymap")
 
 
+# ---------
+# Get GPS time
+# ---------
+runfile = open('bayeswave.run', 'r')
+for line in runfile:
+    spl = line.split()
+    if ' '.join(spl[0:2]) == 'trigger time':
+        trigtime = float(spl[3])
+        print trigtime
+
+gps = LIGOTimeGPS(trigtime)
+sidtime = XLALGreenwichSiderealTime(gps, 0)
+sidtime = sidtime % (np.pi*2)
+
 # -- Try to add a star
-# mdc = ft.Mdc()
-theta, phi = mdc.get_theta_phi(961576872.107)
-# pixel = hp.ang2pix(NSIDE, theta, phi)
+# theta, phi = (2.58350751, -2.782675)
+mdc = ft.Mdc()
+theta, phi = mdc.get_theta_phi(trigtime)
+print theta, phi
+phi += sidtime
+pixel = hp.ang2pix(NSIDE, theta, phi)
 hp.projscatter(theta, phi, s=200, c='yellow', marker=(5,1))
 plt.savefig('skymap.png')
 
@@ -44,7 +64,7 @@ plt.savefig('ra_vs_sample.png')
 
 # -- Attempt to autocorrelate ra
 
-if False:
+if True:
     autocorr = np.correlate(ralist, np.append(ralist, ralist), mode='valid')
     plt.figure()
     plt.plot(autocorr/(autocorr.size))