#!/usr/bin/python
#Standard python libs
import sys
import os
# import datetime
import numpy as np
import h5py 
import matplotlib.pyplot as plt
from math import *
#Libs related to scipy and matplotlib
from scipy import *
from scipy.fftpack import fft
from scipy.fftpack.helper import fftfreq

sys.path.append("./" )
# time_integrator_analysis was created by Jose.
# Jose's function was used directly here.
from time_integrator_analysis import findpeaks, measure_damping, hht_damping_and_shift
from mycolor_fun import *

h5in_file_name=sys.argv[1]
gamma=sys.argv[2]
gamma=float(gamma)
h5file_in=h5py.File(h5in_file_name,"r")

# h5file_in=h5py.File("veri_newmark_dynamic.h5.feioutput","r")

disp=h5file_in['/Model/Nodes/Generalized_Displacements'][()]
time=h5file_in['/time'][()]



# The required displacement.
node_displ=disp[6][:]
last=len(node_displ)-1
len_time=len(time)
node_displ=np.delete(node_displ,[last],None)


dt=time[1]-time[0]

peak_indices=findpeaks(node_displ)

measured_period=time[peak_indices[2]]-time[peak_indices[1]]


alpha=0.0
N = node_displ.shape[0]
D = fft(node_displ[:])
f = fftfreq(N, dt)
xi, fs, Ys = measure_damping(f[0:N/2], abs(D[0:N/2]))	

# **********************************************
# The measured xi is not accurate by the previous method
# Use the other method to measure xi:
peak_indices=findpeaks(node_displ)
pks = []
for index in peak_indices:
	pks.append(node_displ[index])
deltaSize=len(pks)-1;
delta=np.zeros([deltaSize,1]);

for i in xrange(len(delta)):
	delta[i] =1.0/(i+1.0)*np.log(pks[0]/pks[i+1])
deltaAverage=np.mean(delta);
xi_measure=deltaAverage/2./pi;
# print "old xi = " + str(xi)
# print "new xi = " + str(xi_measure)
# **********************************************

peak_index = 0
max_Y = 0
T_system=1.0
for i in arange(len(fs)):
    if (abs(Ys[i]) > max_Y) :