##########################################################################
# program: arbit_rk4.py
# author: Tom Irvine
# version: 1.6
# date: December 15, 2011
# description:  Calculate the response of an SDOF system to a base input.
#               using the Runge-Kutta fourth order method.
#
#               The input file must have two columns: time(sec) & accel(G)
##########################################################################


from tompy import read_two_columns,signal_stats
from tompy import WriteData2,enter_float
from tompy import enter_damping,enter_fn
from tompy import time_history_plot

from numpy import size
from numpy import zeros,linspace,std
from sys import stdin
from math import sqrt

import matplotlib.pyplot as plt

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

def mean_calc(nums):
    return float( sum(nums) ) / len(nums)

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


print "The file must have two columns: time(sec) & accel(G)"

a,b,num =read_two_columns()
b_old=b

t1=a

time_history_plot(a,b_old,1,'Time(sec)','Accel(G)','Base Input','time_history')

sr,dt,mean,sd,rms,skew,kurtosis,dur=signal_stats(a,b,num)

##  sr,dt=sample_rate_check(a,b,num,sr,dt)


iu,fn,omegan,period=enter_fn()


if(iu==1):
    scale=386
else:
    scale=9.81

for i in range(0,num):
        b[i]*=scale



fna=fn

damp,Q=enter_damping()



omegad=omegan*sqrt(1-pow(damp,2))
domegan=damp*omegan

two_domegan=2*domegan

omn2=(omegan**2)

#*******************************************************************************

nt=size(a)-1

#*******************************************************************************

num_fn=1


x=zeros(num)
y=zeros(num)
accel=zeros(num)


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


z0=enter_float()
x[0]=z0


for i in range (0,nt):

    h=t1[i+1]-t1[i]

    hd2=h/2

    a1=           b[i]
    a2= (b[i+1] + b[i])/2.
    a3= a2
    a4=  b[i+1]

    T1=t1[i]
    X1= x[i]
    Y1= y[i]
    F1=-a1-two_domegan*Y1-omn2*X1

    T2=t1[i]+hd2
    X2= x[i]+Y1*hd2
    Y2= y[i]+F1*hd2
    F2=-a2-two_domegan*Y2-omn2*X2

    T3=t1[i]+hd2
    X3= x[i]+Y2*hd2
    Y3= y[i]+F2*hd2
    F3=-a3-two_domegan*Y3-omn2*X3

    T4=t1[i]+h
    X4= x[i]+Y3*h
    Y4= y[i]+F3*h
    F4=-a4-two_domegan*Y4-omn2*X4

    x[i+1]=x[i]+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 )   # relative displacement
    y[i+1]=y[i]+(h/6.)*( F1 + 2*F2 + 2*F3 + F4 )   # relative velocity

    accel[i+1]=-two_domegan*Y4-omn2*X4


if(iu==1):
    acc_resp=accel/386.
else:
    acc_resp=accel/9.81


rd_resp=x

x_pos = max(acc_resp)
x_neg = abs(min(acc_resp))

sss=std(acc_resp)

crest=max(x_pos/sss,x_neg/sss)

mean_acc=mean_calc(acc_resp)

print "\n Acceleration Response (G)"
print "\n  max = %8.4g   min = %8.4g  "      % (x_pos,x_neg)
print "   mean = %8.4g   std dev = %8.4g " % (mean_acc,sss)
print "   crest factor = %8.4g \n"         % (crest)
#


rd_pos = max(rd_resp)
rd_neg = abs(min(rd_resp))
mean_rd=mean_calc(rd_resp)

print "\n Relative Displacement Response (inch) "
print "\n  max = %8.4g   min = %8.4g  "      % (rd_pos,rd_neg)
print "   mean = %8.4g   std dev = %8.4g \n" % (mean_rd,std(rd_resp))
#
#*******************************************************************************
#
#   Velocity - to be added in future revision
#
#*******************************************************************************

iacc=1
ird=1

if(iacc==1):
    print " "
    print "Enter the output acceleration filename: "
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(num,t1,acc_resp,file_path)


if(ird==1):
    print " "
    print "Enter the output relative displacement filename: "
    output_file_path =stdin.readline()
    file_path = output_file_path.rstrip('\n')
    WriteData2(num,t1,rd_resp,file_path)

#*******************************************************************************

print " "
print " view plots"

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

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

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


#*******************************************************************************

if(iacc==1):
    title_string='Acceleration Response fn='+str(fna)+' Hz  Q='+str(Q)
    for i in range(1,200):
        if(Q==float(i)):
            title_string='Acceleration Response fn='+str(fna)+' Hz  Q='+str(i)
            break;

    time_history_plot(t1,acc_resp,2,'Time(sec)','Accel(G)',title_string,'accel_resp')

#*******************************************************************************

if(ird==1):
    title_string='Relative Displacement Response fn='+str(fna)+' Hz  Q='+str(Q)
    for i in range(1,200):
        if(Q==float(i)):
            title_string='Relative Displacement Response fn='+str(fna)+' Hz  Q='+str(i)
            break;

    time_history_plot(t1,rd_resp,3,'Time(sec)','Rel Disp (in)',title_string,'rd_resp')

plt.show()