#! /usr/bin/python

import numpy as np
import math
import matplotlib.pyplot as plt
from  units import *
from param_plot import *

def nextpow2(i):
  n = 2
  while n < i: n = n * 2
  return n
  
#data = np.loadtxt('tracking_expol/2D_probe_velocity_p00001.res')
data = np.loadtxt('tracking_dnt/2D_probe_density_p00000.res')
#print dt*10, 1.0/(dt*10)

time = data[:,0]*dt
p1 = data[:,1]
p2 = data[:,2]
deltaP = (p1-p2)*press_lp
nData = len(time)
vel = data[:,2]
#nSample = 1930
nSample = 2048
shift = 40
#print nSample
nEnd = nData - shift
nStart = nData - nSample - shift


t1 = 0
#t2 = time[nEnd-1]-time[nStart-1]
#freq = nSample/2.0*np.linspace(t1,t2,nSample/2)
#print freq
#interval = time[nEnd-1]-time[nEnd-2]
#print interval
#nSample = len(x[nStart:nEnd])
res = np.fft.fft(vel[nStart:nEnd],nSample)/nSample
#res = np.fft.fft(p1[nStart:nEnd],nSample)/nSample
freq = np.array(range(nSample/2))/(nSample/2.0)*(1.0/2)/interval#/2.0
#freq = np.linspace(0,1,nSample/2)/dt/100.0
#freq = np.fft.fftfreq(nSample,d=interval)
#print dt*10, time[nEnd-1]-time[nEnd-2], data[:,0][nEnd-1]-data[:,0][nEnd-2]
#print freq
#print deltaP[nStart:nEnd]
#print time[nStart:nEnd]
#print res
#print 2.0*np.absolute(res)
#plt.plot(res.real, res.imag,'x')

#print len(freq)
#plt.plot(f,2.0*np.absolute(res[0:nSample/2+1]))
plt.plot(time[nStart:nEnd],deltaP[nStart:nEnd])
#plt.plot(time,deltaP)
res2 = 2.0*np.absolute(res)
#plt.plot(freq[1:50], res2[1:50])
plt.show()

#Fs = 64
#T = 1.0/Fs
#L = 64
#t = np.array(range(0,L))*T
#x = 0.5+0.7*np.sin(2.0*math.pi*5*t) #+ np.sin(2*math.pi*120*t)
#y = x #+ np.random.randn(len(t))
##print t
##print y
#NFFT = nextpow2(L)
##print NFFT
#Y = np.fft.fft(y,NFFT)/L
##Y = fft(y,NFFT)/L
##print Y
##print len(Y)
#f = Fs/2.0*np.linspace(0,1,NFFT/2+1) 
##print NFFT, NFFT/2
##print T
##print f
#ft = np.fft.fftshift(np.fft.fftfreq(NFFT/2, d=T)/2) + NFFT/4
#ft2 = np.fft.fftshift(np.fft.fftfreq(NFFT/2, d=T))
#print ft
#ft2 = np.array(range(NFFT/2))/(NFFT/2.0)*(1.0/2)
#print ft2
#plt.plot(Fs*t,y)
#plt.plot(f,2.*np.absolute(Y[0:NFFT/2+1]))
#plt.plot(ft2,2.0*np.absolute(Y[0:NFFT/2]))
#plt.show()
