##########################################################################
# program: sdof_initial.py
# author: Tom Irvine
# version: 1.4
# date: September 20, 2013
# description:  single-degree-of-freedom system subjected to initial
#               condition excitation via odeint
##########################################################################

from __future__ import print_function
import sys

if sys.version_info[0] == 2:
    print ("Python 2.x")
    import Tkinter as tk
    from tkFileDialog import asksaveasfilename
    import tkMessageBox
           
if sys.version_info[0] == 3:
    print ("Python 3.x")    
    import tkinter as tk 
    from tkinter.filedialog import asksaveasfilename 
    import tkinter.messagebox as tkMessageBox     
    

from sys import stdin

from scipy.integrate import odeint

from tompy import enter_damping,enter_fn,enter_initial
from tompy import time_history_plot
from tompy import GetInteger2,WriteData2

import numpy as np

import matplotlib.pyplot as plt

##########################################################################

def calc_deri(yvec, time, nuc, omc):
   return (yvec[1], -nuc * yvec[1] - omc * yvec[0])

##########################################################################

iu,fn,omegan,period=enter_fn()

damp,Q=enter_damping()

v0,d0=enter_initial(iu)

##########################################################################

omegan2=omegan**2

dur=20*period
dt=period/50
NT=int(round(dur/dt))


cdm = 2*damp*omegan

print (' ')
print (' omegan=%8.4g rad/sec fn=%8.4g Hz' %(omegan,fn))

print (' cdm=%8.4g (rad/sec) ' %(cdm))

time_vec = np.linspace(0, dur, NT)

yarr = odeint(calc_deri, (d0, v0), time_vec, args=(cdm, omegan2))

d=yarr[:,0]
v=yarr[:,1]

a=-cdm*v -omegan2*d

if(iu==1):
    a=a/386
else:
    a=a/9.81

if(fn <=10):
    fna=round(fn,2)

if(fn>10 and fn <=100):
    fna=round(fn,1)

if(fn>100):
    fna=round(fn)

dtitle_string='Displacement Response fn='+str(fna)+' Hz  Q='+str(Q)
vtitle_string='Velocity Response fn='+str(fna)+' Hz  Q='+str(Q)
atitle_string='Acceleration Response fn='+str(fna)+' Hz  Q='+str(Q)

for i in range(1,200):
    if(Q==float(i)):
        dtitle_string='Displacement Response fn='+str(fna)+' Hz  Q='+str(i)
        vtitle_string='Velocity Response fn='+str(fna)+' Hz  Q='+str(i)
        atitle_string='Acceleration Response fn='+str(fna)+' Hz  Q='+str(i)
        break;


if(iu==1):
    time_history_plot(time_vec,d,1,'Time(sec)','Disp(in)',dtitle_string,'disp')
    time_history_plot(time_vec,v,2,'Time(sec)','Vel(in/sec)',vtitle_string,'vel')
else:
    time_history_plot(time_vec,d,1,'Time(sec)','Disp(n)',dtitle_string,'disp')
    time_history_plot(time_vec,v,2,'Time(sec)','Vel(m/sec)',vtitle_string,'vel')

time_history_plot(time_vec,a,3,'Time(sec)','Accel(G)',atitle_string,'accel')


print (" ")
print (" Write output data to files? ")
print (" 1=yes 2=no")

iwd=GetInteger2()

if(iwd==1):

    print (" ")
    print ("Enter the output displacement filename: ")
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(NT,time_vec,d,file_path)

    print (" ")
    print ("Enter the output velocity filename: ")
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(NT,time_vec,v,file_path)


    print (" ")
    print ("Enter the output acceleration filename: ")
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(NT,time_vec,a,file_path)

#########

print (" ")
print (" View plots")
plt.show()