# -*- coding: utf-8 -*-

"""Noise Control Assignment 2
   Jack Bryant, Rylan Norcross, Chuen Ken Ng"""

# Import modules that are needed at a later date

import numpy as np
import matplotlib.pyplot as plt

# Known Parameters
                                    #Units
"""Properties of Air"""
c0 = 343                            #Metres per Second
rho0 = 1.2                          #Kilograms per Metre Cubed

"""Properties of the floor panel (in this case the aluminium plate)"""

width = 2.7                         #Metres
length = 2.5                        #Metres
area = width*length                 #Metres Squared
dampingLF = 0.1                     #Arbitary
thickness = 4e-3                    #Metres
density = 2700.                     #Kilograms per Metre Cubed
massFloor = density*area*thickness  #Kilograms
muFloor = massFloor/thickness       #Kilograms per Metre Squared
youngsModulus = 69e9

"""Properties of the engine"""

massEngine = 1700.                  #Kilograms
resonanceFreqEngine = 6.            #Hertz

"""Experimental data"""

#Frequencies used in Hz
freq = [31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000]
#Angular Frequencies used in rad/s
angularFreq = [2*np.pi*x for x in freq]
#Measured SPL of Engine at each Frequence in dB
engineSPL = [95, 104, 109.5, 110.5, 109, 111, 109.5, 104.5, 98]
#Measured vertical vibration spectrum of engine at each Frequency in dB
engineAcceleration = [143, 146, 147, 146.5, 150.5, 158, 162.5, 163.5, 163]
#Measured noise spectrum due to rolling noise at each Frequency in dB
noiseSPL =[90, 85, 72, 61, 57, 52, 50, 44, 38]

# Sound Reduction Index

def stiffnessPerUnitAreaCalculator(E, L):
    s = E/L
    return s

def soundReductionIndexPlotter(mu, frequency, angularFrequency, rho, c, s, lossfactor):
    angularResonantFrequency = np.sqrt(s/mu)
    resonantFrequency = angularResonantFrequency/(2*np.pi)
    print ""
    print "Resonance Frequency = " + str(resonantFrequency)
    print ""
    soundReductionIndices = []
    for w in angularFrequency:
        fluidLoadingParameter = w * mu / (rho * c)
        if w < angularResonantFrequency:
            soundReductionIndex = 20*np.log10(s/(2*w*rho*c))
            soundReductionIndices.append(soundReductionIndex)
        else:
            if w == angularResonantFrequency:
                # Stiffness controlled
                a = (s * lossfactor / w) / (2 * rho * c)
                b = fluidLoadingParameter
                c = s / (w * rho * c)
                transmissionCoefficient = 1 / ((1 + a) ** 2 + 0.25 * ((b - c) ** 2))
            else:
                transmissionCoefficient = (1 + (fluidLoadingParameter / 2) ** 2) ** -1
            soundReductionIndex = 10*np.log10(1/transmissionCoefficient)
            soundReductionIndices.append(soundReductionIndex)
    plt.figure()
    plt.plot(frequency, soundReductionIndices)
    plt.xscale('log')
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Sound Reduction Index (dB)")
    plt.grid()
    plt.show()

# Rylan

def dBToAccel(AdB=engineAcceleration):
    accelerations = []
    for a in AdB:
        A = (10 ** (a / 10.0)) * 10e-6
        accelerations.append(A)
    return accelerations


def plotAcceleration(frequencies=freq):
    fig1 = plt.figure()
    fig1.add_subplot(111)
    plt.plot(np.log10(frequencies), dBToAccel())
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Acceleration (m/s^2)")
    plt.title("Acceleration")
    plt.grid()
    return None


def accelToVel(radialFreqs=angularFreq, accelerations=dBToAccel()):
    # Assuming A's are peak acceleration values at each frequency f
    velocities = []
    for w, A in zip(radialFreqs, accelerations):
        V = A / w
        velocities.append(V)
    return velocities


def plotVelocity(frequencies=freq):
    fig2 = plt.figure()
    fig2.add_subplot(111)
    plt.plot(np.log10(frequencies), accelToVel())
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Velocity (m/s)")
    plt.title("Veloctiy")
    plt.grid()
    return None

plotAcceleration()
plotVelocity()

def YZbeam(radialFreqs=angularFreq, E=youngsModulus, rho=density):
    b = 75 * 1e-3  # width of beam
    d = 180 * 1e-3  # height of beam
    I = b * (d ** 3.0) / 12
    A = b * d
    mobilities = []
    impedances = []
    for w in radialFreqs:
        zbeam = 2 * (1 + 1j) * (w ** 0.5) * ((E * I) ** (0.25)) * ((rho * A) ** (0.75))
        impedances.append(zbeam)
        ybeam = 1 / zbeam
        mobilities.append(ybeam)
    return mobilities, impedances


def mountsKTotal(fnat=resonanceFreqEngine, m=massEngine):
    """Caclulates the total stiffness of the two mounts"""
    wnat = fnat * 2 * np.pi  # rad/s
    kmounttot = (wnat ** 2.0) * m
    return kmounttot


def YZmounts(radialFreqs=angularFreq):
    # mounts modelled as a spring mobility and Impedance
    mobilities = []
    impedances = []
    k = mountsKTotal()
    for w in radialFreqs:
        Ys = (1j * w) / k
        mobilities.append(Ys)
        Zs = 1 / Ys
        impedances.append(Zs)
    return mobilities, impedances


def YZengine(radialFreqs=angularFreq, m=massEngine):
    # engine modelled as a mass mobility and Impedance
    mobilities = []
    impedances = []
    for w in radialFreqs:
        Ye = -1j / (w * m)
        mobilities.append(Ye)
        Ze = 1 / Ye
        impedances.append(Ze)
    return mobilities, impedances

# Vibration Isolation

def freeVelocity(fnat=resonanceFreqEngine, m=massEngine, radialFreqs=angularFreq, acceleration=engineAcceleration):
    Y,Z = YZengine(radialFreqs, m)
    force = [m*a for a in acceleration]
    freeVelocity = [abs(f*mob) for f,mob in zip(Y,force)]
    return freeVelocity

freeVelocities = freeVelocity()
print ""
print "Free Velocities"
print freeVelocities
print ""

def vibrationPower(fnat=resonanceFreqEngine, m=massEngine, radialFreqs=angularFreq, acceleration=engineAcceleration):
    freeVel = freeVelocity()
    squareFreeVel = [x**2 for x in freeVel]
    rmsFreeVel = (sum(squareFreeVel)/len(squareFreeVel))**0.5
    Yr,Zr = YZbeam(radialFreqs, m)
    W = [np.real(z)*(rmsFreeVel**2) for z in Zr]
    return W

def meanSquareVel(fnat=resonanceFreqEngine, m=massEngine, radialFreqs=angularFreq, acceleration=engineAcceleration, areaLength=10):
    S = 10*width
    power = vibrationPower()
    meanSquareVelocity = [w/(dampingLF*]

# Radiation Efficiency

def radiationEfficiency(rho, c, S, V0, vnormal):
    v = 1/2 * (V0**2)
    W = rho*c*(v)*S
    efficiency = W/(rho*c*S*(vnormal**2))
    return efficiency

# Plot Def

def plots():
    fig1 = plt.figure()
    fig1.add_subplot(111)
    plt.title("Mobility plots for the mounts, beams and engine")
    plt.plot(freq, 20 * np.log10(np.abs(YZbeam()[0])), label='Beam mobility')
    plt.plot(freq, 20 * np.log10(np.abs(YZmounts()[0])), label='Mount mobility')
    plt.plot(freq, 20 * np.log10(np.abs(YZengine()[0])), label='Engine mobility')
    plt.grid()
    plt.legend(fontsize='small', loc='best')
    plt.xlabel('Log Frequency (Hz)')
    plt.ylabel('Mobility (dB)')
    plt.semilogx()
    fig2 = plt.figure()
    fig2.add_subplot(111)
    plt.title("Impedance plots for the mounts, beams and engine")
    plt.plot(freq, 20 * np.log10(np.abs(YZbeam()[1])), label='Beam impedance')
    plt.plot(freq, 20 * np.log10(np.abs(YZmounts()[1])), label='Mount impedance')
    plt.plot(freq, 20 * np.log10(np.abs(YZengine()[1])), label='Engine impedance')
    plt.grid()
    plt.legend(fontsize='small', loc='best')
    plt.xlabel('Log Frequency (Hz)')
    plt.ylabel('Impedance (dB)')
    plt.semilogx()

    plt.show()


# Function Calling

x,y = YZbeam()
print x

stiff = stiffnessPerUnitAreaCalculator(youngsModulus, length)
soundReductionIndexPlotter(muFloor, freq, angularFreq, rho0, c0, stiff, dampingLF)

plots()