import matplotlib.pyplot as plt
import numpy as np
import math
import sys
import xlwt
from scipy import signal

signal_id = [] 
x = [[0 for col in range(5)] for row in range(10000)]
grav = 9806.65
rtod = 180 / np.pi
dtor = np.pi / 180

def FileRead(f):
    
    i = 0
    while True:
        
        line = f.readline()
        if not line:
            f.close()
            break
        
        if i == 0:
            signal_id = line.split()  
            sg_cnt = len(signal_id)
        else: 
            x[i-1][:] = map(float,line.split())
            
        i = i + 1
    
    ndata = i - 1
    print('$ ndata, sg_cnt : ',ndata,sg_cnt)
    print('$ signal_id : ',signal_id)

    return ndata, sg_cnt


def impulse_response(ndata, x, filename):
    
    print('$ impulse_response') 
    
    for i in range(0,ndata):                # Check test data
        if(np.mod(i,100) == 0):
            print('x :',x[i][0],x[i][1],x[i][2],x[i][3],x[i][4])
            

#   Please insert the calculation porcess to obtain gain and phase of Bode plot...
#   Input  Time Signal :  x[:][0],x[:][1] 
#   Output Time Signal :  x[:][0],x[:][2] 

    print(' ')
    print(' ')
    print(' ')
    print(' *** Please help me !! ^^ ')
    print(' ')

#             
            
# In the past, I did fft-convert two time series data of'x0(time):x1' and'x0(time):x2' with FORTRAN code.
# Using these frequency domain spectrum, I calculated the gain and phase of the 'x2 output' for the 'x1 input' and drew a bode plot.

# There are many Bode Plot related functions in SciPy, 
# but I would like to know how to obtain the gain and phase of Bode Plot simply 
# by inputting two time series data obtained from experiments like me!!            
                        
#....................................................................
#            
# scipy.signal.bode
# scipy.signal.bode(system, w=None, n=100)[source]
# Calculate Bode magnitude and phase data of a continuous-time system.

# Parameters
# systeman instance of the LTI class or a tuple describing the system.
# The following gives the number of elements in the tuple and the interpretation:

# 1 (instance of lti)

# 2 (num, den)

# 3 (zeros, poles, gain)

# 4 (A, B, C, D)

# warray_like, optional
# Array of frequencies (in rad/s). Magnitude and phase data is calculated for every value in this array. If not given a reasonable set will be calculated.

# nint, optional
# Number of frequency points to compute if w is not given. The n frequencies are logarithmically spaced in an interval chosen to include the influence of the poles and zeros of the system.

# Returns
# w1D ndarray
# Frequency array [rad/s]

# mag1D ndarray
# Magnitude array [dB]

# phase1D ndarray
# Phase array [deg]            
            
    
#-------------------------------------------------------------------

while True:

    print(' ')
    print('=====================================')
    print('$ Input File Name(*.txt) ?')
    print('=====================================')
    filename = input()
    
    if ("_imp"in filename):
        f = open(filename+".txt", 'r');
        ndata, sg_cnt = FileRead(f)             
        impulse_response(ndata, x, filename)   
                        
    else:
        print('*** Normal End of Job !')
        sys.exit()