# -*- coding: utf-8 -*-
# 2009-12-05 Moritz

import numpy as np
import matplotlib.pyplot as plt
import obspy
from obspy.signal import bandpass
plt.close('all')


tr = obspy.read("data/GR.FUR..LHN.D.2003.268")[0]
t = tr.stats.starttime
tr.trim(t+4600,t+7500) #zoom window
df, npts = tr.stats.sampling_rate, tr.stats.npts
time = np.arange(0,npts)/float(df)


plt.figure(0)
plt.plot(time, tr.data)
plt.title("Original Seismogram")
plt.ylabel("Amplitude")
plt.xlabel("Time [s]")

plt.figure(1)
dat = bandpass(tr.data, 0.25, 0.5, df=df, corners=2)
plt.plot(time, dat)
plt.title('Bandpass T0=%.1fs T1=%.1fs'% (1.0/.25,1.0/.5))
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')


for i in xrange(1,10):
    plt.figure(i+1)
    f0 = 1.0/(10*(i+1))
    f1 = 1.0/(10*i)
    dat = bandpass(tr.data, f0, f1, df=df, corners=2)
    plt.plot(time, dat)
    plt.title('Bandpass T0=%.1fs T1=%.1fs'% (1.0/f0,1.0/f1))
    plt.xlabel('Time [s]')
    plt.ylabel('Amplitude')


plt.show()