#!/usr/bin/env python
from __future__ import division
from gwpy.timeseries import TimeSeries
import numpy
from scipy.signal import medfilt

start_time = 1114644300+14*60
end_time = 1114656395-60

#start_time = 1114585962
#end_time = 1114605616
#end_time = 1114619829-120

cache='snsw.lcf'

wind=TimeSeries.fetch('H1:PEM-EX_WIND_ROOF_WEATHER_MPH.mean,m-trend',start_time,end_time)
st = int(wind.epoch.gps)
et = st+60*len(wind)

range=TimeSeries.read(cache,'H1:DMT-SNSW_EFFECTIVE_RANGE_MPC.max',st,et)

wind_arr = [x.value for x in wind]
range_arr = [x.value for x in range]
#Clean up the dta a bit with a median filter
wind_arr = medfilt(wind_arr,3)
range_arr = medfilt(range_arr,3)

a,b = numpy.polyfit(wind_arr,range_arr,1)
print a,b

windupper = int(max(wind_arr)+3)
min_range,max_range = 7., 10.

from matplotlib import use
use('Agg')
import pylab
pylab.figure()
pylab.title('Wind-range correlation from %u, color=minutes after start' % start_time)
pylab.scatter(wind_arr,range_arr, c=numpy.arange(len(wind_arr)))
xarr=numpy.arange(windupper)
pylab.plot(xarr,a*xarr+b,'k-')
pylab.xlim(0.,windupper)
pylab.ylim(min_range,max_range)
pylab.xlabel('Wind speed (MPH)')
pylab.ylabel('H1 sensemon range (Mpc)')
pylab.colorbar()
pylab.savefig('/home/lundgren/public_html/aDQ/H1-minirun/wind_range.png')

pylab.figure()
pylab.subplot(211)
pylab.plot(numpy.arange(len(wind_arr)), wind_arr)
pylab.ylabel('Wind speed (MPH)')
pylab.subplot(212)
pylab.plot(numpy.arange(len(range_arr)), range_arr)
pylab.ylabel('H1 sensemon range (Mpc)')
pylab.xlabel('Minutes after %u' % start_time)
pylab.savefig('/home/lundgren/public_html/aDQ/H1-minirun/wind_time.png')