from lsl.common import stations

lwa1 = stations.lwa1
stands = lwa1.getStands()

print stands[:3]
print len(stands)

for i in xrange(5):
    print str(stands[i])

antennas = lwa1.getAntennas()

print antennas[:3]
print len(antennas)

for i in xrange(5):
    print str(antennas[i])

for ant in antennas:
    if ant.stand.id == 31 and ant.pol == 0:
        print str(ant.stand)
        print ant.stand.x, ant.stand.y, ant.stand.z

import math

ant31 = None
ant45 = None
for ant in antennas:
    if ant.stand.id == 31 and ant.pol == 0:
        ant31 = ant
        print str(ant.stand)
    elif ant.stand.id == 45 and ant.pol == 0:
        ant45 = ant
        print str(ant.stand)
        
dx,dy,dz = ant31.stand-ant45.stand
print math.sqrt( dx**2 + dy**2 + dz**2 )

for ant in antennas:
    if ant.stand.id == 31:
        print str(ant.stand), ant.pol
        print ant.stand.x, ant.stand.y, ant.stand.z

%matplotlib inline
import numpy
from matplotlib import pyplot as plt

ant = antennas[0]

print "Cable delay for '%s' @ 49 MHz: %.1f ns" % (ant.cable.id, ant.cable.delay(49e6)*1e9,)
print "Cable gain @ 49 MHz: %.4f" % ant.cable.gain(49e6)

freqs = numpy.linspace(10e6, 88e6, 101)
delays = ant.cable.delay(freqs)
gains = ant.cable.gain(freqs)

fig = plt.figure()
axD = fig.add_subplot(1, 2, 1)
axD.plot(freqs/1e6, delays*1e9)
axD.set_xlabel("Frequency [MHz]")
axD.set_ylabel("Delay [ns]")

axG = fig.add_subplot(1, 2, 2)
axG.plot(freqs/1e6, gains)
axG.set_xlabel("Frequency [MHz]")
axG.set_ylabel("Gain")
plt.show()

from lsl.correlator import uvUtils

uv = uvUtils.computeUVTrack(antennas[:32:2])
print uv.shape

fig = plt.figure()
ax = fig.gca()
for i in xrange(uv.shape[0]):
    ax.plot(uv[i,0,:], uv[i,1,:])
ax.set_xlabel('u [$\\lambda$]')
ax.set_ylabel('v [$\\lambda$]')
plt.show()

import time
from lsl import astro
from lsl.misc import geodesy

# Get EOPs for five days ago
jd = astro.unix_to_utcjd(time.time() - 3600*24*5)
mjd = jd - astro.MJD_OFFSET
eop = geodesy.getEOP(mjd)

# Report
print "MJD: %i" % eop.mjd
print "-> x: %.3f\"" % eop.x
print "-> y: %.3f\"" % eop.y
print "-> UT1-UTC Diff.: %.3f s" % eop.utDiff
print "-> Type: %s" % eop.type

# Get EOPs for the last five days
jd1 = astro.unix_to_utcjd(time.time() - 3600*24*5)
jd2 = jd1 + 5
mjd1 = jd1 - astro.MJD_OFFSET
mjd2 = jd2 - astro.MJD_OFFSET
mjds = range(int(mjd1), int(mjd2)+1)
eops = geodesy.getEOP(mjds)

for eop in eops:
    print "MJD: %i" % eop.mjd
    print "-> x: %.3f\"" % eop.x
    print "-> y: %.3f\"" % eop.y
    print "-> UT1-UTC Diff.: %.3f s" % eop.utDiff
    print "-> Type: %s" % eop.type