Entry 4085

from scipy.io import matlab as mio
from pylab import plot, show, ylabel, xlabel, title, figure, legend, annotate
import numpy as np

def enforce_complex_sym(array, nyquist):
	'''Enforce complex symmetry'''
	array[:,nyquist+1:] = np.conj(array[:,nyquist-1:0:-1])

def analyse_time_data(t, a, plot_title='', variable='', plot_label='', yunits='m/s^{2}'):
	plot(t, a, label=plot_label)
	if len(plot_title)>0: title(plot_title)
	ylabel('%s $\ (%s)$' % (variable, yunits))
	xlabel('Time (s)')
	max_coord = (t[a.argmax()], a[a.argmax()])
	annotate('max %s' % variable, max_coord, (max_coord[0] *1.5, max_coord[1]*1.1), 
				arrowprops=dict(arrowstyle="simple",
                                fc="0.6", ec="none",
                                connectionstyle="arc3,rad=0.3")   )

def plot_responses(t, arrayN, var, unit, N=3):
	figure()
	for i in range(N):
		analyse_time_data(t, arrayN[i,:], variable=var, plot_label='Floor %d' % (1+i), yunits=unit)
	legend()
	title('%s of each floor in Kobe earthquake' % var)


# Load earthquake data from matlab file
kobe_data = mio.loadmat('Kobe.mat')
a, t, dt = [kobe_data[index][0] for index in ['f', 't', 'dt']]

analyse_time_data(t, a, 'Kobe Earthquake Data', 'Ground Acceleration')

# data on each story:
m = 10000
k = 1600000
c = 13000

# Construct System matricies:
M = np.matrix(np.diag(3*[m]))

C = np.matrix( [[2*c,  -c,  0],
				[ -c, 2*c, -c],
				[  0,  -c,  c]])
					
K = np.matrix( [[2*k, -k,  0],
				[-k, 2*k, -k],
				[0,   -k,  k]])

F = ( -M*np.matrix(np.ones(3)).transpose() ) * a
N = F.shape[1]
nyquist = (N/2)

# Transform into frequency domain
Fs = np.fft.fft(F, axis=1)
w = np.array(2 * np.pi * np.fft.fftfreq(N, dt))

# Calculate and apply the transfer function upto nyquist frequency
Vs, dVs, ddVs = [np.zeros([3,N],np.complex) for i in range(3)]
for i in xrange(nyquist):
	Gs = ((K - (w[i])**2 *M) + 1j*w[i]*C).I
	Vs[:,i] = (Gs * np.c_[Fs[:,i]]).reshape(3)
	dVs[:,i] =  Vs[:,i] * 1j * w[i]
	ddVs[:,i] = Vs[:,i] * w[i]**2

[enforce_complex_sym(array, nyquist) for array in [Vs, dVs, ddVs]]
# Convert back to the time domain
vt, dvt, ddvt = [np.fft.ifft(array, axis=1).real for array in [Vs, dVs, ddVs]]

# Plot the response
plot_responses(t, vt, 'Displacement', 'm')
plot_responses(t, dvt, 'Velocity', 'm/s')
plot_responses(t, ddvt, 'Acceleration', 'm/s^{2}')

# Show all the plots
show()