################################################################################
# program: arbit_force.py
# author: Tom Irvine
# Email: tom@vibrationdata.com
# version: 1.1
# date: September 12, 2013
# description:  Calculate the response of an SDOF system to an applied force.
#               The file must have two columns: time(sec) & force
################################################################################

from __future__ import print_function

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

from scipy.signal import lfilter
from numpy import zeros,std

from math import pi,exp,sqrt,cos,sin

import matplotlib.pyplot as plt

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

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 natural frequency (Hz)")
            
fn=enter_float()            

omegan=2*pi*fn

stiffness=mass*omegan**2

print (" ")
if(iu==1):
    print (" stiffness=%8.4g lbf/in" %stiffness)
else:
    print (" stiffness=%8.4g N/m" %stiffness)    
print (" ")


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)
sample_rate_check(a,b,num,sr,dt)

print (sr)

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

ac=zeros(3,'f')
bc=zeros(3,'f') 

omegad=omegan*sqrt(1.-(damp**2))
domegan=damp*omegan
 
cosd=cos(omegad*dt)
sind=sin(omegad*dt)

domegadt=domegan*dt

eee1=exp(-domegadt)
eee2=exp(-2.*domegadt)

ecosd=eee1*cosd
esind=eee1*sind 

omeganT=omegan*dt

phi=(2*damp)**2-1

a1= 2.*ecosd
a2=-eee2


ac[0]=1
ac[1]=-a1
ac[2]=-a2

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

    
DD1=(omegan/omegad)*phi
DD2=2*DD1
    
df1=2*damp*(ecosd-1) +DD1*esind +omeganT
df2=-2*omeganT*ecosd +2*damp*(1-eee2) -DD2*esind
df3=(2*damp+omeganT)*eee2 +(DD1*esind-2*damp*ecosd)

MD=(mass*omegan**3*dt)
df1=df1/MD
df2=df2/MD
df3=df3/MD


bc[0]=df1
bc[1]=df2
bc[2]=df3


disp_resp=lfilter(bc, ac, b, axis=-1, zi=None)
#
x_pos = max(disp_resp)
x_neg = abs(min(disp_resp))
    
sss=std(disp_resp)
    
crest=max(x_pos/sss,x_neg/sss)    
    
mean_disp=mean_calc(disp_resp)  
    
print ("\n Displacement Response (in)")
print ("\n   max = %8.4g   min = %8.4g  "      % (x_pos,x_neg))
print ("   mean = %8.4g   std dev = %8.4g " % (mean_disp,sss))

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

VV1=damp*omegan/omegad

vf1=(-ecosd+VV1*esind)+1
vf2=eee2-2*VV1*esind-1
vf3=ecosd+VV1*esind-eee2

VD=(mass*omegan**2*dt)
vf1=vf1/VD
vf2=vf2/VD
vf3=vf3/VD

bc[0]=vf1
bc[1]=vf2
bc[2]=vf3

vel_resp=lfilter(bc, ac, b, axis=-1, zi=None)
#
x_pos = max(vel_resp)
x_neg = abs(min(vel_resp))
    
sss=std(vel_resp)
    
crest=max(x_pos/sss,x_neg/sss)    
    
mean_vel=mean_calc(vel_resp)  
    
print ("\n Velocity Response (in/sec)")
print ("\n   max = %8.4g   min = %8.4g  "      % (x_pos,x_neg))
print ("   mean = %8.4g   std dev = %8.4g " % (mean_vel,sss))

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

af1=esind/(mass*omegad*dt)
af2=-2*af1
af3=af1

bc[0]=af1
bc[1]=af2
bc[2]=af3

acc_resp=lfilter(bc, ac, b, axis=-1, zi=None)
#

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

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 (" ")
print (" view plots")

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

title_string='Applied Force'

if(iu==1):  
    time_history_plot(t,b,1,'Time(sec)','Force(lbf)',title_string,'force')
else:
    time_history_plot(t,b,1,'Time(sec)','Force(N)',title_string,'force')
    

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

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

if(iu==1):  
    time_history_plot(t,disp_resp,2,'Time(sec)','Disp(in)',title_string,'d_resp')
else:
    time_history_plot(t,disp_resp,2,'Time(sec)','Disp(m)',title_string,'d_resp')
    
#*******************************************************************************

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

if(iu==1):    
    time_history_plot(t,vel_resp,3,'Time(sec)','Vel(in/sec)',title_string,'vel_resp')
else:
    time_history_plot(t,vel_resp,3,'Time(sec)','Vel(m/sec)',title_string,'vel_resp')
    
#*******************************************************************************

title_string='Acceleration Response fn='+str(fn)+' Hz  Q='+str(Q)
for i in range(1,200):
    if(Q==float(i)):
        title_string='Acceleration Response fn='+str(fn)+' Hz  Q='+str(i)
        break  
    
time_history_plot(t,acc_resp,4,'Time(sec)','Accel(G)',title_string,'accel_resp')


plt.show()