##########################################################################
# program: arbit_force_rk4.py
# author: Tom Irvine
# version: 1.0
# date: December 7, 2011
# description:  Calculate the response of an SDOF system to an applied
#               force using the Runge-Kutta fourth order method.
#
#               The input file must have two columns: time(sec) & force
##########################################################################


from tompy import read_two_columns,signal_stats
from tompy import GetInteger2,enter_float
from tompy import enter_damping
from tompy import time_history_plot

from numpy import size
from numpy import zeros
from math import sqrt,pi

import matplotlib.pyplot as plt

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


print " "
print " Select units "
print " 1=English  2=metric"
iu=GetInteger2()

if(iu==1):
    print(" Enter mass (lbm)")
else:
    print(" Enter mass (kg)")

mass=enter_float()

if(iu==1):
    mass/=386

if(iu==1):
    print(" Enter stiffness (lbf/in)")
else:
    print(" Enter stiffness (N/m)")
            
stiffness=enter_float()            

omegan=sqrt(stiffness/mass)
fn=omegan/(2*pi)

period=1/fn        

damp,Q=enter_damping()


print " "
print " fn=%8.4g Hz" %fn
print " "

if(iu==1):
    print "The file must have two columns: time(sec) & force(lbf)"
else:
    print "The file must have two columns: time(sec) & force(N)"
            

a,b,num =read_two_columns()

t=a
f=b

b_old=b

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

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


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)	

for i in range (1,nt):
    f[i]/=mass

for i in range (1,nt):
    
    h=t[i+1]-t[i]

    hd2=h/2

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

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

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

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

    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


d=x
v=y

a=-two_domegan*v-omn2*d+f

time_vec=t

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


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 "\n Acceleration Response (G)"
print "   max = %8.4g   min = %8.4g  "      % (max(a),min(a))

print " "
print " view plots"

plt.show()