# This may take a bit...
import os
import urllib2
from lsl.reader import tbw

if not os.path.exists('temp.tbw'):
    fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/firstlight/firstLightTBW.dat')
    fh2 = open('temp.tbw', 'wb')
    fh2.write(fh1.read(tbw.FrameSize*1000))
    fh1.close()
    fh2.close()
    
print "TBW Size: %.1f kB" % (os.path.getsize('temp.tbw')/1024.)

from lsl.reader import tbw, errors

fh = open('temp.tbw', 'rb')
for i in xrange(5):
    try:
        frame = tbw.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    print "DP Stand Number: %i" % frame.parseID()
    print "-> Frame Count: %i" % frame.header.frameCount
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f" % frame.data.xy.mean()
fh.close()

from lsl.common.stations import lwa1
from lsl.reader import tbw, errors

# Get the antenna list
antennas = lwa1.getAntennas()

# Read in the first 5 frames
fh = open('temp.tbw', 'rb')
for i in xrange(5):
    try:
        frame = tbw.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    # Calculate the digitizer
    dig = 2*(frame.parseID()-1) + 1
    
    print "DP Stand Number: %i" % frame.parseID()
    print "-> LWA1 Stand: %i" % antennas[dig-1].stand.id
    print "-> Frame Count: %i" % frame.header.frameCount
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f" % frame.data.xy.mean()
fh.close()

# This may take a bit...
import os
import urllib2
from lsl.reader import tbn

if not os.path.exists('temp.tbn'):
    fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/TBN/055987_000496599_DEFAULT_TBN_00000_30.dat')
    fh2 = open('temp.tbn', 'wb')
    fh2.write(fh1.read(tbn.FrameSize*1000))
    fh1.close()
    fh2.close()
    
print "TBN Size: %.1f kB" % (os.path.getsize('temp.tbn')/1024.)

from lsl.reader import tbn, errors

# Read in the first 5 frames
fh = open('temp.tbn', 'rb')
for i in xrange(5):
    try:
        frame = tbn.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    print "DP Stand Number: %i, pol. %i" % frame.parseID()
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f + %.3f j" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)
fh.close()

from lsl.common.stations import lwa1
from lsl.reader import tbn, errors

# Download a small bit of TBN data (1000 frames to be exact)
# This may take a bit...
import os
import urllib2
if not os.path.exists('temp.tbn'):
    fh1 = urllib2.urlopen('http://fornax.phys.unm.edu/lwa/data/TBN/055987_000496599_DEFAULT_TBN_00000_30.dat')
    fh2 = open('temp.tbn', 'wb')
    fh2.write(fh1.read(tbn.FrameSize*1000))
    fh1.close()
    fh2.close()
    
# Get the antenna list
antennas = lwa1.getAntennas()

# Read in the first 5 frames
fh = open('temp.tbn', 'rb')
for i in xrange(5):
    try:
        frame = tbn.readFrame(fh)
    except errors.syncError:
        continue
    except errors.eofError:
        break
        
    # Calculate the digitizer
    stand,pol = frame.parseID()
    dig = 2*(stand-1) + pol + 1
    
    print "DP Stand Number: %i, pol. %i" % (stand,pol)
    print "-> LWA1 Stand, pol.: %i, %i" % (antennas[dig-1].stand.id, antennas[dig-1].pol)
    print "-> Time tag: %.6f s" % frame.getTime()
    print "-> Mean value: %.3f + %.3f j" % (frame.data.iq.mean().real, frame.data.iq.mean().imag)
fh.close()

%matplotlib inline
from lsl.misc.mathutil import to_dB
from matplotlib import pyplot as plt

fh = open('temp.tbn', 'rb')
srate = tbn.getSampleRate(fh)   # Data sample rate in Hz
frame = tbn.readFrame(fh)
cFreq = frame.getCentralFreq()  # Data tuning in Hz
frame.data.iq.shape = (1,frame.data.iq.size)  # SpecMaster expects 2-D data

from lsl.correlator import fx as fxc
freq, spec = fxc.SpecMaster(frame.data.iq, LFFT=512, SampleRate=srate, CentralFreq=cFreq) 
print freq.shape, spec.shape

fig = plt.figure()
ax = fig.gca()
ax.plot((freq-freq.mean())/1e3, to_dB(spec[0,:]))
ax.set_xlabel('Relative Freq [kHz]')
ax.set_ylabel('PSD [arb. dB]')
plt.show()

import math
import ephem
from lsl.common import stations

lwa1 = stations.lwa1
lwaObserver = lwa1.getObserver(2455548.38787, JD=True)
jove = ephem.Jupiter()
jove.compute(lwaObserver)
print "Jupiter:  az -> %.1f, el -> %.1f" % (jove.az*180/math.pi, jove.alt*180/math.pi)

lwaObserver.date = '2010/12/17 21:18:34'

cyga = ephem.FixedBody()
cyga._ra = '19:59:28.30'
cyga._dec = '+40:44:02'
cyga.compute(lwaObserver)
print "Cygnus A:  az -> %.1f, el -> %.1f" % (cyga.az*180/math.pi, cyga.alt*180/math.pi)

from lsl.misc import beamformer

antennas = [lwa1.getAntennas()[0],]

beamdata = beamformer.phaseAndSum(antennas, frame.data.iq, sampleRate=1e5, azimuth=10.9, elevation=83.2)
print beamdata.shape