###########################################################################
# program: srs_conv.py
# author: Tom Irvine
# version: 1.0
# date: October 10, 2011
# description:  Calculate the shock response spectrum for an SDOF system
#               using the convolution intergral method.
#
#               The file must have two columns: time(sec) & accel(G)
###########################################################################

from tompy import read_two_columns,signal_stats,sample_rate_check
from tompy import GetInteger2,GetInteger3,WriteData3
from tompy import enter_damping
from tompy import enter_float
from tompy import time_history_plot
from tompy import srs_plot_pn

from numpy import zeros,concatenate
from numpy import exp,sqrt,cos,sin
from numpy import array,size

from sys import stdin
from math import pi

import matplotlib.pyplot as plt


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

t,b,num =read_two_columns()

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

sr,dt=sample_rate_check(t,b,num,sr,dt)

print(" ")
damp,Q= enter_damping()

print(" ")
print(" Enter starting frequency (Hz) ")        
f1 = enter_float()

print(" ")
print(" Select octave spacing")
print("   1= one-third")
print("   2= one-sixth")
print("   3= one-twelfth")

noct=GetInteger3()
        
if(noct==1):
    oct=1./3.
if(noct==2):
    oct=1./6.
if(noct==3):
    oct=1./12.
    

print(" ")        
print("  include residual:  1=yes 2=no ")

ires = GetInteger2()
        
    
print(" ")
                
print(" ")
print(" Select output metrics")
print(" ")
print("  acceleration:  1=yes 2=no ")

iacc = GetInteger2()
        
print(" ")        
print("  pseudo velocity:  1=yes 2=no ")

ipv = GetInteger2()

print(" ")        
print("  relative displacement:  1=yes 2=no ")

ird = GetInteger2()

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

nt=int(0.8*(1/f1)/dt)

b_old=b
if(ires==1):
    tz=zeros(nt)
    d=concatenate((b,tz))
    b=d
#*******************************************************************************

resp=zeros(num)
fn=zeros(1000)
ac=zeros(3)
bc=zeros(3)


fn[0]=f1    
for j in range(1,999):
    fn[j]=fn[j-1]*(2.**oct)
    if  fn[j] > sr/8.:
        break

num_fn=j-1

x_pos=zeros(num_fn)
x_neg=zeros(num_fn)
pv_pos=zeros(num_fn)
pv_neg=zeros(num_fn)
rd_pos=zeros(num_fn)
rd_neg=zeros(num_fn)


hrd=zeros(num,'f')
gamma=zeros(num,'f')


t=array(t)
t=t-t[0]

if(iacc==1):

    print 'num_fn=%d' % num_fn

    print " "
    print "   fn      Positive    Negative"
    print "  (Hz)     Accel(G)    Accel(G) "
    
    for j in range(0,num_fn):
   
        omegan=2.*pi*fn[j]
        omegad=omegan*sqrt(1.-(damp**2))
        cosd=cos(omegad*dt)
        sind=sin(omegad*dt)
        domegadt=damp*omegan*dt
        
        arg=omegad*t
#
        alpha=(omegad**2-(damp*omegan)**2)/omegad;
        beta =2*damp*omegan        
#
        domegan=damp*omegan

        S= alpha*sin(arg) + beta*cos(arg)        

      
      
        ha=zeros(num,'f')
        ha= exp(-domegan*t)*S

#
        resp=zeros(num,'f');
        w=zeros(num,'f');   

#
#   Convolution with single loop
#
        for i in range (0,num):
            w[0:num-i]=b[i]*ha[0:num-i]
            resp[i:num]+=w[0:num-i]
	    
        resp*=dt;

#
        x_pos[j]= max(resp)
        x_neg[j]= abs(min(resp))
        ns=size(x_pos)

        print "%6.4g   %6.4g   %6.4g " %(fn[j],x_pos[j],x_neg[j])
        
if(ipv==1 or ird==1):       
    
    for j in range(0,num_fn):
   
        omegan=2.*pi*fn[j]
        omegad=omegan*sqrt(1.-(damp**2))
        cosd=cos(omegad*dt)
        sind=sin(omegad*dt)
        domegadt=damp*omegan*dt
        
        arg=omegad*t
#
        alpha=(omegad**2-(damp*omegan)**2)/omegad;
        beta =2*damp*omegan        
#
        domegan=damp*omegan

        S= sin(arg)   
      
        hrd=zeros(num,'f')
        hrd= -exp(-domegan*t)*S/omegad

#
        resp=zeros(num,'f');
        w=zeros(num,'f');   

#
#   Convolution with single loop
#
        for i in range (0,num):
            w[0:num-i]=b[i]*hrd[0:num-i]
            resp[i:num]+=w[0:num-i]
	    
        resp*=dt;

#
        rd_pos[j]= max(resp)*386
        rd_neg[j]= abs(min(resp))*386
        
        pv_pos[j]= omegan*rd_pos[j]
        pv_neg[j]= omegan*rd_neg[j]
        
#
#*******************************************************************************
#
print " "
print " begin plots "
print " "

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

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

if(iacc==1):
    
    plt.figure(2)
    ns=size(x_pos)
    fp=fn[0:ns]
    srs_plot_pn(1,1,fp,x_pos,x_neg,damp,'accel_srs_plot')  
    

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

if(ipv==1):
    plt.figure(3)
    plt.plot(fp, pv_pos, label="positive")
    plt.plot(fp, pv_neg, label="negative")
    plt.xscale('log')
    plt.yscale('log')
    plt.grid(True)
#
    title_string='Pseudo Velocity Shock Response Spectrum Q='+str(Q)
#
    for i in range(1,200):
        if(Q==float(i)):
            title_string='Pseudo Velocity  Shock Response Spectrum Q='+str(i)
            break;
#    
    plt.title(title_string)
    plt.xlabel('Natural Frequency (Hz) ')
    plt.ylabel('Peak Velocity (in/sec)')
    plt.grid(True, which="both")
    plt.savefig('srs_plot')
    plt.legend(loc="upper left")

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

if(ird==1):
    plt.figure(4)
    plt.plot(fp, rd_pos, label="positive")
    plt.plot(fp, rd_neg, label="negative")
    plt.xscale('log')
    plt.yscale('log')
    plt.grid(True)
#
    title_string='Relative Displacement Shock Response Spectrum Q='+str(Q)
#
    for i in range(1,200):
        if(Q==float(i)):
            title_string='Relative Displacement Shock Response Spectrum Q='+str(i)
            break;
#    
    plt.title(title_string)
    plt.xlabel('Natural Frequency (Hz) ')
    plt.ylabel('Peak Rel Disp (in)')
    plt.grid(True, which="both")
    plt.savefig('srs_plot')
    plt.legend(loc="upper left")

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

plt.show()