##########################################################################
# program: sdof_base_sine.py
# author: Tom Irvine
# version: 1.5
# date: January 13, 2012
# description:  single-degree-of-freedom system subjected to
#               sinusoidal base excitation
##########################################################################

from sys import stdin

from scipy.integrate import odeint

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

import numpy as np

from math import pi
import matplotlib.pyplot as plt

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

def calc_deri(yvec, time, nuc, omc, amp, omega_base):
   f=-nuc*yvec[1] -omc*yvec[0] -amp*np.sin(omega_base*time)
   return (yvec[1],f)

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

# iu=units   1=English  2=metric

iu,fn,omegan,period=enter_fn()

damp,Q=enter_damping()



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


cdm = 2*damp*omegan
omegan2=omegan**2

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

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

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

print ' '
print(" Enter base acceleration amplitude (G)")

amp=enter_float()


if(iu==1):
    amp*=386.
else:
    amp*=9.81



print " "
if(iu==1):
    print " Enter initial relative displacement(inch) "
else:
    print " Enter initial relative displacement(m) "

d0=enter_float()

v0=0


print(" Enter base excitation frequency (Hz)")

fbase=enter_float()

omega_base=fbase*(2*pi)

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

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

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

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

a=np.zeros(NT,float)
base=np.zeros(NT,float)

for i in range (0,NT):
    a[i]=-cdm*v[i] -omegan2*d[i]
    base[i]=amp*np.sin(omega_base*time_vec[i])

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

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

iwd=GetInteger2()

if(iwd==1):

    print " "
    print "Enter the output relative 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 acceleration filename: "
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(NT,time_vec,a,file_path)

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

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='Relative Displacement Response fn='+str(fna)+' Hz  Q='+str(Q)
vtitle_string='Relative Velocity Response fn='+str(fna)+' Hz  Q='+str(Q)
atitle_string='Absolute Acceleration Response fn='+str(fna)+' Hz  Q='+str(Q)

for i in range(1,200):
    if(Q==float(i)):
        dtitle_string='Relative Displacement Response fn='+str(fna)+' Hz  Q='+str(i)
        vtitle_string='Relative Velocity Response fn='+str(fna)+' Hz  Q='+str(i)
        atitle_string='Absolute 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(m)',dtitle_string,'disp')
    time_history_plot(time_vec,v,2,'Time(sec)','Vel(m/sec)',vtitle_string,'vel')

plt.figure(3)
plt.plot(time_vec, a, label="response")
plt.plot(time_vec, base, label="input")
plt.xlabel('Time (sec)')
plt.ylabel('Accel (G)')
plt.grid(True)
plt.title(atitle_string)
plt.savefig('accel')
plt.legend(loc="upper right") 
plt.draw()
    

Amax=max(a)
Amin=min(a)

RDmax=max(d)
RDmin=min(d)


print " "
print " Acceleration Response (G) "
print "   max=%8.4g    min=%8.4g \n" %(Amax,Amin)
print " "

if(iu==1):
    print " Relative Displacement Response (in) "
else:
    print " Relative Displacement Response (m) "

print "   max=%8.4g    min=%8.4g  \n" %(RDmax,RDmin)    
print " "

print " view plots"

plt.show()