Final

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 30 13:10:14 2018

@author: Lawrence
"""
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('Agg')
from obspy import read, Trace
from obspy.signal import PPSD
from obspy.io.xseed import Parser
from obspy.signal.invsim import corn_freq_2_paz
import time
from obspy.signal.invsim import corn_freq_2_paz
from obspy.io.sac.sacpz import attach_resp
from obspy.imaging.cm import pqlx
import numpy as np
import math
import matplotlib.mlab as mlab

paz_sts2 = {
    'poles': [-3.691000E-02 + 3.702000E-02j, -3.691000E-02 - 3.702000E-02j, -3.430000E+02 + 0j, -3.700000E+02 + 4.670000E+02j,
              -3.700000E+02 - 4.670000E+02j, -8.360000E+02 + 1.522000E+03j, -8.360000E+02 - 1.522000E+03j, -4.900000E+03 + 4.700000E+03j, -4.900000E+03 - 4.700000E+03j,-6.900000E+03 + 0j, -1.500000E+04 + 0j],
    'zeros': [0, 0, 0, -3.920000E+02 + 0j, -1.960000E+03 + 0j, -1.490000E+03 + 1.740000E+03j, -1.490000E+03 - -1.740000E+03j],
    'gain':  4.344928E+17 * 800000,
    'sensitivity': 7.543000E+02}
paz_1hz=corn_freq_2_paz(1.0, damp=0.707)
paz_1hz['sensitivity']= 1.0
#Inputs:
sta = 'GNGN'    #Station name
input_fils_A= "C:/Users/lawre/Desktop/Thesis_Research/New_Data/Centaur-6_1270_20180429_190000.miniseed" #Make sure that "ls" will return files
input_fils_B= "C:/Users/lawre/Desktop/Thesis_Research/New_Data/Centaur-6_1281_20180429_190000.miniseed" #Second Seismometer
win_len = 3600. #Duration of each segment in seconds ppsd_length
ovr_lap = 0.5  #Overlap is the overlap between consecutive segments
sens_resp = "C:/Users/lawre/Desktop/Thesis_Research/S-ConfigurationSensitivityFilesforexistingTrilliumSensor/BHZ.TrilliumCompact.120.754.dataless" #Sensor Response
digi_resp = "C:/Users/lawre/Desktop/Thesis_Research/S-ConfigurationSensitivityFilesforexistingTrilliumSensor/XX.NN443..HHZ.CENTAUR.2.250.OFF.dataless"#Digitizer Response
#=======================================================================================
sA = read(input_fils_A)
for chn in ['HHZ']: #Loop through channels
    scur = sA.select(channel=chn) #Select the channel
print sA
resample_rate=10
sA.resample(resample_rate)
print sA
sA_orig = sA.copy() #Copy of original data
sA.simulate(paz_remove=paz_sts2, paz_simulate=paz_1hz)
#sA_orig.plot()
sA[0].plot()

sB = read(input_fils_B)
for chn in ['HHZ']: #Loop through channels
    scur = sB.select(channel=chn) #Select the channel
print sB
sB.resample(resample_rate)
print sB
sB_orig = sB.copy() #Copy of original data
sB.simulate(paz_remove=paz_sts2, paz_simulate=paz_1hz)
#sB_orig.plot()
sB[2].plot()

tr_A = np.array(sA[0])#number 2 is HHZ
tr_B = np.array(sB[2])
distance =np.array(3.0)
tr_difference = np.subtract(tr_A, tr_B)
tr_tilt = np.divide(tr_difference,distance)
times=[]
i=0
while i < len(tr_tilt):
    newval=float((i)*(1./resample_rate))
    times.append(newval)
    i+=1
matplotlib.pyplot.plot(times,tr_tilt)
matplotlib.pyplot.show()

fs = resample_rate
n=tr_difference.size
freqall = np.fft.rfftfreq(n, d=1./fs)
freqplot=freqall
power = np.fft.rfft(tr_tilt, norm="ortho")#the normalising fixed the amplitude to match Josh's result. Applying this changes normalisation from 1/n to 1/sqrt(n)
powerplot =np.absolute(power)
print(freqplot.size)
print(powerplot.size)
A = plt.plot([0.01, 0.2, 0.5], [3e-10, 6e-9, 1e-10],'r-', lw=2, zorder=2)
plt.plot(freqplot,powerplot, zorder=1)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Frequency (Hz)')
plt.ylabel('rad/sqrt(Hz)')
plt.show()