chan='L1:LSC-POP_A_RF9_I_ERR_DQ'
st=1162024217
dur=600

data=TimeSeries.fetch(chan,st,st+dur)

srate=data.sample_rate.value

"""
Explanation: Here's some data that sees the caduceus / raygun noise.
End of explanation
"""


twid=int(srate) # Length of short spectra
zpadsec=8 # Length (sec) for oversampled FFT
zpadwid=int(srate*zpadsec)

win1=sig.hann(twid)
def myfft(data,tt):
    idx=int(tt*srate)-twid/2
    tmp=zeros(int(zpadsec*srate))
    tmp[:twid]=win1*data.value[idx:idx+twid]
    tmp=rfft(tmp)
    return tmp[100*zpadsec:140*zpadsec]

def track_lines(data,tt):
    fft1=myfft(data,tt)
    i1=abs(fft1[:20*zpadsec+1]).argmax()
    i2=abs(fft1[20*zpadsec:]).argmax()+20*zpadsec
    return (tt,i1/float(zpadsec)+100.,i2/float(zpadsec)+100) #,fft1[i1],fft1[i2])

"""
Explanation: The functions below take a short segment of data, do a zero-padded FFT (to do some extra smoothing), then find the largest peak above and below 120 Hz.
End of explanation
"""


tst=array([track_lines(data,tt) for tt in arange(1,dur-1,0.25)])

"""
Explanation: Take the FFT and find peaks every quarter of a second.
End of explanation
"""


filt=sig.firwin(int(8*srate),[105./nyq,135./nyq],window='hann',pass_zero=False)
filt.resize(len(data))
ffilt=abs(rfft(filt))

fdata=ffilt*rfft(data.detrend().value)
fdata[0]=0.
newdata=TimeSeries(irfft(fdata)[int(4*srate):-int(4*srate)],sample_rate=srate)

p1=newdata[:int(srate*120)].spectrogram2(1,0.98).plot()
plt.plot(tst[:,0]-4,tst[:,1],c='r')
plt.plot(tst[:,0]-4,tst[:,2],c='r')
#p1.gca().set_yscale('log')
p1.gca().set_ylim(100,140)

"""
Explanation: To make the spectrogram look nice, we'll do some bandpassing. This is left over from my previous attempt to do a Hilbert transform and get the phase directly, but I like the results.
End of explanation
"""


data2=TimeSeries.fetch('L1:LSC-SRCL_IN1_DQ',st,st+dur)

srate2=data2.sample_rate.value
nyq2=srate2/2.
filt2=sig.firwin(int(8*srate2),[4./nyq2,30./nyq2],window='hann',pass_zero=False)
filt2.resize(len(data2))
ffilt2=abs(rfft(filt2))

fdata2=ffilt2*rfft(data2.detrend().value)
fdata2[0]=0.
newdata2=TimeSeries(irfft(fdata2[:1+512*int(data2.duration.value)])[4*1024:-4*1024],sample_rate=1024)

sgram2=newdata2.spectrogram2(1,0.95)

p1=sgram2.plot()
plt.plot(tst[:,0]-4,1.15*(tst[:,2]-120),c='r')
#p1.gca().set_yscale('log')
p1.gca().set_ylim(5,15)
p1.gca().set_xlim(0,120)

p1=sgram2.plot()
plt.plot(tst[:,0]-4,1.15*(tst[:,2]-120),c='r')
#p1.gca().set_yscale('log')
p1.gca().set_ylim(5,15)
p1.gca().set_xlim(240,360)

"""
Explanation: Now we look at the scattered light arches in a related degree of freedom, SRCL, to see if they're related to the raygun.
End of explanation
"""


