First attempt at healpix skymap

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

postprocess/skymap/make_sky.sh

+#!/bin/sh
+
+echo "# ra sin_dec" > skymap.txt
+chainFullName="$(ls chains/*signal_1_wavechain.dat.0)"
+awk '{print $2"   "$3}' $chainFullName >> skymap.txt
\ No newline at end of file

postprocess/skymap/skygram.py

+
+import numpy as np
+import healpy as hp
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+
+
+plt.close('all')
+NSIDE = 64
+print "The number of pixels is {0}".format(hp.nside2npix(NSIDE))
+
+# -------------
+# Read sky coordinates
+# -------------
+filename = 'skymap.txt'
+ralist, sin_dec = np.loadtxt('skymap.txt', unpack=True)
+dec = np.arcsin(sin_dec)
+thetalist = np.pi/2.0 - dec
+
+# -- Make a histogram
+m = np.zeros(hp.nside2npix(NSIDE))
+for theta, ra in zip(thetalist, ralist):    
+    pixel =  hp.ang2pix(NSIDE, theta, ra)
+    m[pixel] += 1
+hp.mollview(m,title="Sweet skymap")
+plt.savefig('skymap.png')
+
+
+# -- Look for autocorrelations
+plt.figure()
+plt.plot(ralist)
+plt.xlabel('Sample')
+plt.ylabel('RA')
+plt.savefig('ra_vs_sample.png')
+
+# -- Attempt to autocorrelate ra
+autocorr = np.correlate(ralist, ralist, mode='full')
+auto = autocorr[autocorr.size/2:]
+plt.figure()
+plt.plot(auto)
+
+plt.ylabel('Autocorrelation')
+plt.xlabel('Number of sample offsets')
+plt.savefig('RA_autocorr.png')
+