SEA model examples

Two rooms separated by a wall

# -*- coding: utf-8 -*-
"""
Example script for lecture example
for Ulf Orenius
"""
import numpy as np
import matplotlib.pyplot as plt

# pyva packages
import pyva.models as mds
import pyva.coupling.junctions as con

import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC


import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC

import pyva.useful as uf

mycolors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf', '#17becf', '#17becf']
 
# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels(range='SEA')
# Print option
printsw = True
#printsw = False
figpath = '../source/images/'
plt.close('all')

# wall dimensions
Ly = 4.
Lz = 2.5
S = Lz*Ly

# Additional room dimensions
Lx1 = 3.
Lx2 = 5.
# Absorption area
As1  = 8.
As2  = 10. # As2 = S

# wall thickness
h = 0.05

# junction properties
angle_S = 0
angle_R = 90*np.pi/180

# Create materials
concrete = matC.IsoMat(E=3.8e9,nu=0.33,rho0=1250.)
air = matC.Fluid()

# Create properties
concrete_5cm = stPC.PlateProp(h,concrete)

# Create plate subsystems
wall  = st2Dsys.RectangularPlate(2, Ly,Lz,prop=concrete_5cm, eta = 0.03)
room1 = ac3Dsys.RectangularRoom(1, Lx1, Ly, Lz, air, absorption_area = As1, damping_type= ['surface'] )
# use receiving room with absorption
room2 = ac3Dsys.RectangularRoom(3, Lx2, Ly, Lz, air, absorption_area = As2, damping_type= ['surface'])

J123 = con.AreaJunction((room1,wall,room2))#,area_wall)
#om = 2*np.pi*np.logspace(2,np.log10(10000),3)


# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)
om    = omega.data
freq  = om/2/np.pi 

# define load
# room power
power1mWatt = lC.Load(omega, 0.001*np.ones(omega.shape), dof.DOF(1,0,dof.DOFtype(typestr = 'power')), name = '1Watt')

# wall power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
two_rooms = mds.HybridModel((wall,room1,room2),xdata=omega)
#connect both
two_rooms.add_junction({'areaJ_12':J123})
two_rooms.add_load('1mWatt',power1mWatt) # add 1Watt per band 

#pdof = J123.get_wave_DOF()

# Calculate specific CLF factors
eta_13 = J123.CLF(om,(0,2))
eta_31 = J123.CLF(om,(2,0))
eta_12 = J123.CLF(om,(0,1),)
eta_21 = J123.CLF(om,(1,0))

eta_all = J123.junction_matrix(om)

# check tau
tau_diff = J123.non_res.transmission_diffuse(om,signal = False)



#%% plot the different modal densities
plt.figure(1)
plt.plot(freq,wall.modal_density(om,3),label = 'wall')
plt.plot(freq,room1.modal_density(om),label = 'room1')
plt.plot(freq,room2.modal_density(om),label = 'room2')
# plt.loglog(f,n_1,':',label = 'n1 VA1')
# plt.loglog(f,n_2,':',label = 'n2 VA1')
# plt.loglog(f,n_3,':',label = 'n3 VA1')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$n(\omega)/$Hz$^{-1}$')
plt.legend()
plt.show

if printsw:
     plt.savefig(figpath+'two_rooms_modal_density.png')


#%% plot modal overlap 
plt.figure(2)
plt.plot(freq,wall.modal_overlap(om,3),label = 'wall B')
plt.plot(freq,wall.modal_overlap(om,5),label = 'wall LS')
plt.plot(freq,room1.modal_overlap(om),label = 'room1')
plt.plot(freq,room2.modal_overlap(om),label = 'room2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$M$')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_modal_overlap.png')

#%% plot modes in band 
plt.figure(3)
plt.plot(freq,wall.modes_in_band(om,3),label = 'wall B')
plt.plot(freq,wall.modes_in_band(om,5),label = 'wall LS')
plt.plot(freq,room1.modes_in_band(om),label = 'room1')
plt.plot(freq,room2.modes_in_band(om),label = 'room2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('modes in band')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_modes_in_band.png')


 



#%% PLot radation efficiency
sigma = wall.radiation_efficiency(om,fluid = air)
sigma_simple = wall.radiation_efficiency_simple(om,fluid = air)

# Show coincidence frequency
f_c = concrete_5cm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

plt.figure(4)
plt.loglog(freq,sigma,label = '[Lep1982]')
plt.loglog(freq,sigma_simple,label = 'ISO EN 12354-1')
#plt.loglog(om_VA1,sigma_VA1,label = 'VAOne')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\sigma$')
#plt.ylim(0,0.005)
#plt.xticks(2*np.pi*fc,fclabels)
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_radiation_efficiency.png')


#%% Plot CLFs
plt.close(5)
plt.figure(5)
eta_13.plot(5,yscale = 'log',fulllegstr = ['$\eta_{13}$'],ls='-')
eta_31.plot(5,yscale = 'log',fulllegstr = ['$\eta_{31}$'],cs=[mycolors[1]],ls = '--')
eta_12.plot(5,yscale = 'log',fulllegstr = ['$\eta_{1,2B}$'],cs=[mycolors[2]],ls = ':',marker = [' '])
eta_21.plot(5,yscale = 'log',fulllegstr = ['$\eta_{2B,1}$'],cs=[mycolors[3]],ls = '-.',marker = [' '])
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\eta$')
plt.ylim(1E-8,1E-2)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])
plt.xscale('log')

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_etas.png')

#%% Solve model
two_rooms.create_SEA_matrix(sym = 1)
two_rooms.solve()

# Derive paths of power input
pow_in_room2 = two_rooms.power_input(3)

#%% Energry result plots
two_rooms.energy.plot(20,xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('wall B','wall LS','room1','room2'))
plt.xlabel('$f_c/$Hz')
plt.ylabel('$E/$ J')
#plt.ylim(bottom=0.1)
plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_energy.png')



#%% Engineering units velocity
plt.figure(30,figsize=(5,5/4*3))
two_rooms.result.plot(30,ID=[2],dof=[3],xscale = 'log',yscale = 'log',
                       fulllegstr = ('wall B',))
#plt.plot(om_VA1,v_wall_B,':',label = 'VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$v_{\\rm rms}/$m s$^{-1}$')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_velocity.png')

#%% engineering units pressure
plt.figure(31,figsize=(5,5/4*3))

two_rooms.result.plot(31,ID=[1,3],xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('room 1','room 2'))
#plt.plot(om_VA1,p_room1,':',label = 'room1 VA1')
#plt.plot(om_VA1,p_room2,':',label = 'room2 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_pressure.png')
    
# Calculate transmission loss from pressure difference
# Note: This is allowed because wall and absorption surface are equal S=A2
p1 = two_rooms.result[2].ydata.flatten()
p2 = two_rooms.result[3].ydata.flatten()
tau = (p2/p1)**2

#%% Plot transmission loss 

plt.figure(32)
plt.semilogx(freq,-10*np.log10(tau),label = 'wall TL')
#plt.semilogx(f,20*np.log10(tau_VA1),':',label = 'pure VA1')
#plt.semilogx(freq,-10*np.log10(mass_law(om)),':',label = 'mass law')
#plt.semilogx(freq,20*np.log10(wall.prop.mass_per_area*freq)-48,':',label = 'mass law diffuse')
plt.semilogx(freq,-10*np.log10(tau_diff),':',label = 'mass law TMM')
plt.xticks(fc[2:],fclabels[2:])
plt.xlabel('$f_c/$Hz')
plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_TL.png')

#%% Plot input power
    
pow_in_room2.plot(33,yscale = 'log',fulllegstr = ['$\Pi_{room1}$','$\Pi_{wall}$'],xscale = 'log')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_power_in.png')






    

Two rooms with floor separated by a wall

# -*- coding: utf-8 -*-
"""
Example script for lecture example
for Ulf Orenius
"""
import numpy as np
import matplotlib.pyplot as plt

# pyva packages
import pyva.models as mds
import pyva.coupling.junctions as con

import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC


import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC

import pyva.useful as uf

mycolors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf', '#17becf', '#17becf']
 
# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels(range='hybrid')
# Print option
printsw = True
printsw = False
figpath = '../source/images/'

plt.close('all')

# wall dimensions
Ly = 4.
Lz = 2.5
S = Lz*Ly

# Additional room dimensions
Lx1 = 3.
Lx2 = 5.
# Absorption area
As1  = 8.
As2  = 10. # As2 = S

# wall thickness
h   = 0.05
h_f = 0.17 

# junction properties
angle_S = 0
angle_R = 90*np.pi/180

# Create materials
concrete = matC.IsoMat(E=3.8e9,nu=0.33,rho0=1250.)
air = matC.Fluid()

# Create properties
concrete_5cm   = stPC.PlateProp(h,concrete)
concrete_17cm  = stPC.PlateProp(h_f,concrete)

# Create plate subsystems
wall   = st2Dsys.RectangularPlate(2, Ly,Lz,prop=concrete_5cm, eta = 0.03)
floor1 = st2Dsys.RectangularPlate(4, Lx1,Ly,prop=concrete_17cm, eta = 0.03)
floor2 = st2Dsys.RectangularPlate(5, Lx2,Ly,prop=concrete_17cm, eta = 0.03)
room1 = ac3Dsys.RectangularRoom(1, Lx1, Ly, Lz, air, absorption_area = As1, damping_type= ['surface'] )
# use receiving room with absorption
room2 = ac3Dsys.RectangularRoom(3, Lx2, Ly, Lz, air, absorption_area = As2, damping_type= ['surface'])

# Area Junctions
J123 = con.AreaJunction((room1,wall,room2))
J14  = con.AreaJunction((room1,floor1))
J35  = con.AreaJunction((room2,floor2))

# Line Junctions
J425 = con.LineJunction((floor1,wall,floor2),length = Ly, thetas = (0,90,180))


# Frequency range
omega = mC.DataAxis.octave_band(f_min=2*np.pi*50,f_max=2*np.pi*10000)
om    = omega.data
freq  = om/2/np.pi 

# define load
# room power
# power1mWatt = lC.Load(omega, 0.001*np.ones(omega.shape), dof.DOF(1,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
force10Nrms = lC.Load(omega, 10*np.ones(omega.shape), dof.DOF(4,3,dof.DOFtype(typestr = 'force')), name = '10N')

# wall power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
two_rooms = mds.HybridModel((wall,room1,room2,floor1,floor2),xdata=omega)
#connect all
two_rooms.add_junction({'areaJ_123':J123})
two_rooms.add_junction({'areaJ_14':J14})
two_rooms.add_junction({'areaJ_35':J35})
two_rooms.add_junction({'lineJ_425':J425})

#two_rooms.add_load('10Watt',power1mWatt) # add 1Watt per band 
two_rooms.add_load('10N',force10Nrms)# add force excitatio to wave_DOF 3 of system 4
#
# Calculate specific CLF factors
eta_13 = J123.CLF(om,(0,2))
eta_31 = J123.CLF(om,(2,0))
eta_12 = J123.CLF(om,(0,1),)
eta_21 = J123.CLF(om,(1,0))

eta_all = J123.junction_matrix(om)

# check tau
tau_diff = J123.non_res.transmission_diffuse(om,signal = False)

# read reference results (keep if ulf will check with VA1)
#f,n_1,n_2,n_3 = np.loadtxt ('2rooms_m_dens.csv',
#                    unpack = True,
#                    usecols = (0,2,1,5), skiprows = 1,
#                    delimiter = ',')

# f,sigma_VA1 = np.loadtxt ('2rooms_sigma.csv',
#                     unpack = True,
#                     usecols = (0,1), skiprows = 1,
#                     delimiter = ',')

# f,p_room1,p_room2 = np.loadtxt ('2rooms_p_rms.csv',
#                     unpack = True,
#                     usecols = (0,1,2), skiprows = 1,
#                     delimiter = ',')

# f,v_wall_B = np.loadtxt ('2rooms_v_rms.csv',
#                     unpack = True,
#                     usecols = (0,1), skiprows = 1,
#                     delimiter = ',')

# f,pow_wall,pow_room1 = np.loadtxt ('2rooms_pow_in.csv',
#                     unpack = True,
#                     usecols = (0,3,2), skiprows = 1,
#                     delimiter = ',')

# f,tau_VA1,tau_eff_VA1 = np.loadtxt ('2rooms_tau.csv',
#                     unpack = True,
#                     usecols = (0,1,2), skiprows = 1,
#                     delimiter = ',')

# f,eta_13_VA1,eta_31_VA1,eta_12_VA1,eta_21_VA1,eta_23_VA1,eta_32_VA1 = np.loadtxt ('2rooms_eta.csv',
#                     unpack = True,
#                     usecols = (0,1,2,6,5,3,4), skiprows = 1,
#                     delimiter = ',')


# om_VA1 = f*2*np.pi


#%% plot the different modal densities
plt.figure(1)
plt.plot(freq,wall.modal_density(om,3),label = 'wall')
plt.plot(freq,floor1.modal_density(om,3),label = 'floor1')
plt.plot(freq,floor2.modal_density(om,3),label = 'floor2')
# plt.loglog(f,n_1,':',label = 'n1 VA1')
# plt.loglog(f,n_2,':',label = 'n2 VA1')
# plt.loglog(f,n_3,':',label = 'n3 VA1')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$n(\omega)/$Hz$^{-1}$')
plt.legend()
plt.show

if printsw:
     plt.savefig(figpath+'two_rooms_floor_modal_density.png')


#%% plot modal overlap 
plt.figure(2)
plt.plot(freq,wall.modal_overlap(om,3),label = 'wall B')
#plt.plot(freq,wall.modal_overlap(om,5),label = 'wall LS')
plt.plot(freq,floor1.modal_overlap(om,3),label = 'floor1')
plt.plot(freq,floor2.modal_overlap(om,3),label = 'floor2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$M$')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_modal_overlap.png')

#%% plot modes in band 
plt.figure(3)
plt.plot(freq,wall.modes_in_band(om,3),label = 'wall B')
#plt.plot(freq,wall.modes_in_band(om,5),label = 'wall LS')
plt.plot(freq,floor1.modes_in_band(om,3),label = 'floor1')
plt.plot(freq,floor2.modes_in_band(om,3),label = 'floor2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('modes in band')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_modes_in_band.png')


 



#%% PLot radation efficiency
sigma1 = floor1.radiation_efficiency(om,fluid = air)
sigma2 = floor2.radiation_efficiency(om,fluid = air)

# Show coincidence frequency
f_cw = concrete_5cm.coincidence_frequency(air.c0)/2/np.pi
f_cf = concrete_17cm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency wall = {0}'.format(f_cw))
print('coincidence frequency floor = {0}'.format(f_cf))

plt.figure(4)
plt.semilogx(freq,sigma1,label = 'floor1')
plt.semilogx(freq,sigma2,label = 'floor2')
#plt.loglog(om_VA1,sigma_VA1,label = 'VAOne')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\sigma$')
#plt.ylim(0,0.005)
#plt.xticks(2*np.pi*fc,fclabels)
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_radiation_efficiency.png')


#%% plot4
def mass_law(omega,theta=0):
    z0 = air.z0
    m  = wall.prop.mass_per_area
    return 1/(1+(m*omega/2/z0)**2)

def eta_f(omega,V,area):
    c0 = air.c0
    return c0*area/4/V/omega*mass_law(omega)

plt.figure(5)
eta_13.plot(5,yscale = 'log',fulllegstr = ['$\eta_{13}$'],ls='-')
eta_31.plot(5,yscale = 'log',fulllegstr = ['$\eta_{31}$'],cs=[mycolors[1]],ls = '--')
#plt.loglog(om_VA1,eta_13_VA1,':',label = '13 VA1')
#plt.loglog(om_VA1,eta_31_VA1,':',label = '31 VA1')
# #%% plot5
#eta_12.plot(4,yscale = 'log',fulllegstr = ['$\eta_{1,2B}$'],cs=[mycolors[2]],ls = ':',marker = [' '])
#eta_21.plot(4,yscale = 'log',fulllegstr = ['$\eta_{2B,1}$'],cs=[mycolors[3]],ls = '-.',marker = [' '])
#plt.loglog(om_VA1,eta_12_VA1,':',label = '13 VA1')
#plt.loglog(om_VA1,eta_21_VA1,':',label = '31 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\eta$')
#plt.ylim(0,0.005)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])
plt.xscale('log')

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_etas.png')

#%% Solve model
two_rooms.create_SEA_matrix(sym = 1)
two_rooms.solve()

# Derive paths of power input
pow_in_room1 = two_rooms.power_input(1)
pow_in_room2 = two_rooms.power_input(3)


#%% Enrgry result plots
two_rooms.energy.plot(20,xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                        fulllegstr = ('wall B','wall LS','room1','room2','floor1 B','floor1 LS','floor2 B','floor2 LS'))
plt.xlabel('$f_c/$Hz')
plt.ylabel('$E/$ J')
#plt.ylim(bottom=0.1)
plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_energy.png')



#%% Engineering units velocity
plt.figure(30,figsize=(5,5/4*3))
two_rooms.result.plot(30,ID=[2,4,5],dof=[3,3,3],xscale = 'log',yscale = 'log',
                       fulllegstr = ('wall B','floor1 B','floor2 B'))
#plt.plot(om_VA1,v_wall_B,':',label = 'VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$v_{\\rm rms}/$m s$^{-1}$')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_velocity.png')

#%% engineering units pressure
plt.figure(31,figsize=(5,5/4*3))

two_rooms.result.plot(31,ID=[1,3],xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('room 1','room 2'))
#plt.plot(om_VA1,p_room1,':',label = 'room1 VA1')
#plt.plot(om_VA1,p_room2,':',label = 'room2 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_pressure.png')
    
# Calculate transmission loss from pressure difference
# Note: This is allowed because wall and absorption surface are equal S=A2
p1 = two_rooms.result[2].ydata.flatten()
p2 = two_rooms.result[3].ydata.flatten()
tau = (p2/p1)**2


#%% Plot input power
    
pow_in_room1.plot(32,yscale = 'log',xscale = 'log')
#plt.loglog(om_VA1,pow_wall,':')
#plt.loglog(om_VA1,pow_room1,':')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_power_in1.png')

pow_in_room2.plot(33,yscale = 'log',xscale = 'log')
#plt.loglog(om_VA1,pow_wall,':')
#plt.loglog(om_VA1,pow_room1,':')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_power_in2.png')





    

Box cover of sound source

# -*- coding: utf-8 -*-
"""
Pyva box example without isolation
"""

import numpy as np
import matplotlib.pyplot as plt

import pyva.models as mds
import pyva.coupling.junctions as jun
import pyva.properties.structuralPropertyClasses as stPC
import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC

import pyva.useful as uf

# my packages

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC
import pyva.properties.materialClasses as matC

plt.close('all')

# x-axis tics for better readability
fc,fclabels = uf.get_3rd_oct_axis_labels()
fc,fclabels = fc[1:],fclabels[1:]

# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)

# Plate dimensions
Lx = 1.2
Ly = 1
Lz = 1.1

# Box dimensions
V = Lx*Ly*Lz
A = 2*(Lx*Ly+Ly*Lz+Lx*Lz)
P = 4*(Lx+Ly+Lz)

# Plates thickness
t = 2.0E-3;

# junction angle
angle_R = (0.,90*np.pi/180)

# Create materials
steel = matC.IsoMat(E=210e9,nu=0.3,rho0=7800, eta = 0.0)
air   = matC.Fluid() 

# Create props
steel2mm = stPC.PlateProp(t,steel)


area_dof = dof.DOF(0,0,dof.DOFtype(typestr='area'))


# %% create models
f_c = steel2mm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

# Create plate subsystems
plate1 = st2Dsys.RectangularPlate(1,Lx,Lz,prop = steel2mm)
plate2 = st2Dsys.RectangularPlate(2,Ly,Lz,prop = steel2mm)
plate3 = st2Dsys.RectangularPlate(3,Lx,Lz,prop = steel2mm)
plate4 = st2Dsys.RectangularPlate(4,Ly,Lz,prop = steel2mm)
plate5 = st2Dsys.RectangularPlate(5,Lx,Ly,prop = steel2mm)

room     = ac3Dsys.Acoustic3DSystem(6, V , A, P, air)

# create semi infinite fluids
sif1 = jun.SemiInfiniteFluid((room,plate1), air)
sif2 = jun.SemiInfiniteFluid((room,plate2), air)
sif3 = jun.SemiInfiniteFluid((room,plate3), air)
sif4 = jun.SemiInfiniteFluid((room,plate4), air)
sif5 = jun.SemiInfiniteFluid((room,plate5), air)

juncs =  { 'j61' : jun.AreaJunction((room,plate1)), 
           'j62' : jun.AreaJunction((room,plate2)),
           'j63' : jun.AreaJunction((room,plate3)),
           'j64' : jun.AreaJunction((room,plate4)),
           'j65' : jun.AreaJunction((room,plate5)),
           'j12' : jun.LineJunction((plate1,plate2), Lz, angle_R),
           'j23' : jun.LineJunction((plate2,plate3), Lz, angle_R),
           'j34' : jun.LineJunction((plate3,plate4), Lz, angle_R),
           'j41' : jun.LineJunction((plate4,plate2), Lz, angle_R),
           'j15' : jun.LineJunction((plate1,plate5), Lx, angle_R),
           'j25' : jun.LineJunction((plate2,plate5), Ly, angle_R),
           'j35' : jun.LineJunction((plate3,plate5), Lx, angle_R),
           'j45' : jun.LineJunction((plate4,plate5), Ly, angle_R)
            }



#om = 2*np.pi*np.logspace(2,np.log10(10000),3)




# define load
# room power
power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(6,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
# plate power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
box = mds.HybridModel((plate1,plate2,plate3,plate4,plate5,room),xdata=omega)

#connect both
box.add_junction(juncs)
box.add_SIF({'sif1' : sif1, 
             'sif2' : sif2, 
             'sif3' : sif3, 
             'sif4' : sif4, 
             'sif5' : sif5})

box.add_load('1Watt',power1Watt) # add 1mWatt per band 


#%% solving
box.create_SEA_matrix(sym = 1)
box.solve()

#%% plotting 1

box.result.plot(1,ID = [1,5],xscale = 'log',yscale = 'log',fulllegstr=('a)','a)'))
plt.figure(1)
plt.yscale('log')
plt.ylim(1e-5,0.01)
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
#plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

#%% plotting 2
 
box.result.plot(2,ID = [6],xscale = 'log',yscale = 'linear',fulllegstr=('a)',))   
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
plt.tight_layout()
plt.legend()
plt.show()

# plt.savefig('../source/images/box_pure_pressure.png')

#%% more info
sif1_in = box.power_input('sif1')
sif2_in = box.power_input('sif2')
sif3_in = box.power_input('sif3')
sif4_in = box.power_input('sif4')
sif5_in = box.power_input('sif5')

sif_all = sif1_in.sum()+sif2_in.sum()+sif3_in.sum()+sif4_in.sum()+sif5_in.sum()


#%% plot 3

sif1_in.plot(3,xscale='log',yscale='linear')
plt.figure(3)
#plt.plot(om_VA1,power_in_1_res,':',label = 'VA1 res')
#plt.plot(om_VA1,power_in_1_nonres,':',label = 'VA1 non-res')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_pure_SIF_power_in.png')


#%% plot 4

sif_all.plot(4,xscale='log',res = 'dB',fulllegstr=('a)',))
plt.figure(4)
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
#plt.yticks([40,60,80,100,],[40,60,80,100,])
plt.ylim((60,110))
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_pure_SIF_total.png')

#%% plot 5    
NR = 1/sif_all.ydata[0,:]

plt.figure(5)
plt.semilogx(omega.data,10*np.log10(NR),label = 'a)')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$IL/$dB')
plt.legend()
plt.tight_layout()
plt.show()




    

Box cover with noise control treatment

# -*- coding: utf-8 -*-
"""
Pyva box example with isolation

"""
import numpy as np
import matplotlib.pyplot as plt

import pyva.models as mds
import pyva.coupling.junctions as con
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC

import pyva.useful as uf

# my packages

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC
import pyva.properties.materialClasses as matC

plt.close('all')

# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels()
fc,fclabels = fc[1:],fclabels[1:]

# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)
om    = omega.data

# Plate dimensions
Lx = 1.2
Ly = 1
Lz = 1.1

# Room dimensions
V = Lx*Ly*Lz
A = 2*(Lx*Ly+Ly*Lz+Lx*Lz)
P = 4*(Lx+Ly+Lz)


# Plate thickness
t = 2.0E-3;

# junction properties
angle_R = (0.,90*np.pi/180)

# Create materials
steel = matC.IsoMat(E=210e9,nu=0.3,rho0=7800, eta = 0.0)
air   = matC.Fluid()

# Fibre material
rho_bulk = 0.98*1.20 + 30.
fibre1 = matC.EquivalentFluid(porosity = 0.98, \
                               flow_res = 25000.,\
                               tortuosity = 1.02, \
                               length_visc = 90.e-6, \
                               length_therm = 180.e-6,\
                               rho_bulk = rho_bulk , \
                               rho0 = 1.208, \
                               dynamic_viscosity = 1.81e-5 )

# Thickness of fibre layer
h1 = 0.05

# Create props
steel2mm = stPC.PlateProp(t,steel)
# Create Layer

fibre_5cm = iL.FluidLayer(h1,fibre1)
heavy_1kg = iL.MassLayer(0.001, 1000)
# and noise control treatement
nct_mass = mds.TMmodel((fibre_5cm,heavy_1kg)) 
# change order for absorption
nct_mass4abs = mds.TMmodel((heavy_1kg,fibre_5cm))
nct          = mds.TMmodel((fibre_5cm,)) 



# absorption due to trim
alpha_nct = nct.absorption_diffuse(om,in_fluid=air)
alpha_nct_mass = nct_mass4abs.absorption_diffuse(om,in_fluid=air)

# %% plot abs fig 10 
alpha_nct.plot(10,xscale='log',yscale='linear',fulllegstr = ('5cm fibre',))
alpha_nct_mass.plot(10,xscale='log',yscale='linear',fulllegstr = ('5cm fibre + 1kg',), ls = '--')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\\alpha_s$')
plt.xticks(2*np.pi*fc,fclabels)
plt.tight_layout()
plt.legend()
plt.show()

# %% create absorption
area_dof = dof.DOF(0,0,dof.DOFtype(typestr='area'))
abs_area = alpha_nct*(Lx*Ly)+alpha_nct_mass*(A-Lx*Ly)
# This is a signal those dof must be set to the
abs_area.dof = area_dof

# %% plot 1

abs_area.plot(1,xscale = 'log')

plt.savefig('../source/images/box_isolated_abs_area.png')

# %% plates

f_c = steel2mm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

# Create plate subsystems
plate1 = st2Dsys.RectangularPlate(1,Lx,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate2 = st2Dsys.RectangularPlate(2,Ly,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate3 = st2Dsys.RectangularPlate(3,Lx,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate4 = st2Dsys.RectangularPlate(4,Ly,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate5 = st2Dsys.RectangularPlate(5,Lx,Ly,prop = steel2mm,trim=(nct_mass,'none'))

room   = ac3Dsys.Acoustic3DSystem(6, V , A, P, air,\
                        absorption_area = abs_area ,\
                        damping_type= ['eta','surface']    )



# use receiving room with abso
sif1 = con.SemiInfiniteFluid((room,plate1), air)
sif2 = con.SemiInfiniteFluid((room,plate2), air)
sif3 = con.SemiInfiniteFluid((room,plate3), air)
sif4 = con.SemiInfiniteFluid((room,plate4), air)
sif5 = con.SemiInfiniteFluid((room,plate5), air)


juncs =  { 'j61' : con.AreaJunction((room,plate1)), 
           'j62' : con.AreaJunction((room,plate2)),
           'j63' : con.AreaJunction((room,plate3)),
           'j64' : con.AreaJunction((room,plate4)),
           'j65' : con.AreaJunction((room,plate5)),
           'j12' : con.LineJunction((plate1,plate2), Lz, angle_R),
           'j23' : con.LineJunction((plate2,plate3), Lz, angle_R),
           'j34' : con.LineJunction((plate3,plate4), Lz, angle_R),
           'j41' : con.LineJunction((plate4,plate2), Lz, angle_R),
           'j15' : con.LineJunction((plate1,plate5), Lx, angle_R),
           'j25' : con.LineJunction((plate2,plate5), Ly, angle_R),
           'j35' : con.LineJunction((plate3,plate5), Lx, angle_R),
           'j45' : con.LineJunction((plate4,plate5), Ly, angle_R)
           }

# define load
# room power
power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(6,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
# plate power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')


#create SEA model
box = mds.HybridModel((plate1,plate2,plate3,plate4,plate5,room),xdata=omega)

#connect both
box.add_junction(juncs)
box.add_SIF({'sif1' : sif1, 
             'sif2' : sif2, 
             'sif3' : sif3, 
             'sif4' : sif4, 
             'sif5' : sif5})

box.add_load('1Watt',power1Watt) # add 1Watt per band 



#%% Solve
box.create_SEA_matrix(sym = 1)
box.solve()


#%% plotting 3

box.result.plot(2,ID = [1,5],xscale = 'log',yscale = 'log',fulllegstr=('a)','a)'))
plt.figure(2)
plt.yscale('log')
plt.ylim(1e-5,0.01)
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
#plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

#%% plotting 3
 
box.result.plot(3,ID = [6],xscale = 'log',yscale = 'linear',fulllegstr=('a)',))   

plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
plt.tight_layout()
plt.legend()
plt.show()

#%% more info
sif1_in = box.power_input('sif1')
sif2_in = box.power_input('sif2')
sif3_in = box.power_input('sif3')
sif4_in = box.power_input('sif4')
sif5_in = box.power_input('sif5')

sif_all = sif1_in.sum()+sif2_in.sum()+sif3_in.sum()+sif4_in.sum()+sif5_in.sum()


#%% plot 4

sif1_in.plot(4,xscale='log',yscale='linear')
plt.figure(4)
#plt.plot(om_VA1,power_in_1_res,':',label = 'VA1 res')
#plt.plot(om_VA1,power_in_1_nonres,':',label = 'VA1 non-res')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#%% plot 5

sif_all.plot(5,xscale='log',res = 'dB',fulllegstr=('b)',))
plt.figure(5)
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
#plt.yticks([40,60,80,100,],[40,60,80,100,])
#plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_isolated_SIF_total.png')



#%% plot 6    
NR = 1/sif_all.ydata[0,:]

plt.figure(6)
plt.semilogx(om,10*np.log10(NR),label = 'b)')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$IL/$dB')
plt.legend()
plt.tight_layout()
plt.show()




    

Price and Crocker double panel transmission loss using SEA

# -*- coding: utf-8 -*-
"""
Script related to the paper from Price and Crocker [1] on the sound transmission of double panels.

[1] A. J. Price and M. J. Crocker, “Sound Transmission through Double Panels Using 
    Statistical Energy Analysis, ”The Journal of the Acoustical Society of America, 
    vol. 47, no. 3A, pp. 683–693, Mar. 1970, doi: 10.1121/1.1911951.
"""

import numpy as np
import matplotlib.pyplot as plt

# pyva main model module 
import pyva.models as mds

# Coupling module
import pyva.coupling.junctions as con    # junctions module

# Property modules
import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC

# Layer module for transfermatrix methods
import pyva.systems.infiniteLayers as iL

# SEA system modules
import pyva.systems.structure2Dsystems as st2Dsys # plate systems module
import pyva.systems.acoustic3Dsystems as ac3Dsys  # cavity systems moduke 
import pyva.loads.loadCase as lC                  # loadcase module

# Modules for system matrices and degree of freedom
import pyva.data.dof as dof
import pyva.data.matrixClasses as mC

# Pyva utilities
import pyva.useful as uf

aprop  = dict(facecolor='black', shrink=0.01, width=1 )

#%% Model settings

# Create frequency axis labels
fc,fclabels = uf.get_3rd_oct_axis_labels(range='SEA')
plt.close('all')

#%% Frequency range
omega = mC.DataAxis.octave_band(f_min=2*np.pi*100,f_max=2*np.pi*10000,bands_per_octave = 3)
# Create according ndarrays for plotting
om    = omega.data
freq  = om/2/np.pi 

# Panel dimensions
L1 = 1.55   # panel/DW cavity width
L2 = 0.071  # DW cavity thickness
L3 = 1.97   # panel/DW cavity height

L1_fig3 = 2.2 # height of cavity from fig 3

S = L1*L3   # panel surface
S_dwc_edge = L2*2*(L1+L3)# DW cavity edge surface

# Cavity volumes
V_dwc  = L1*L2*L3 # Volume DW cavity
V1 = 100 # room 1
V2 = 100 # room 2

# Absorption areas in large cavities
As1  = 2. # arbitrary room absorption area
As2  = S  # secong room absorption area equal to junction area for log10(S/As2)=0

# Materials
alu = matC.IsoMat(eta=0) # Alu is the default isomaterial
air = matC.Fluid(eta=0)  # Air is the default fluid

# Plate properties
alu_3_18mm = stPC.PlateProp(0.00318,alu) # 0.85 g/cm^2
alu_6_36mm = stPC.PlateProp(0.00636,alu) # 1.71 g/cm^2
alu_0_88mm = stPC.PlateProp(0.00088,alu) # 0.24 g/cm^2
alu_9_26mm = stPC.PlateProp(0.00925,alu) # 2.5  g/cm^2

#%% Tables for DW cavity damping
omega_min = 2*np.pi*100     # in all arguments pyva uses angular frequency
omega_max = 2*np.pi*10000
freq_tab_BI = mC.DataAxis.octave_band(omega_min,omega_max,typeID = 21)
 
# Create table until 10000Hz
eta_DWC_data = np.zeros(np.shape(freq_tab_BI))
# Fill data until 2000 Hz
eta_DWC_data[0,:14] = 0.01*np.array([[3.18, 2.29, 2.38, 2.89, 2.68, 
                         2.72, 2.79, 3.27, 3.09, 2.80,
                         2.38, 2.00, 1.60, 1.27 ]])

# Fill data from 3200 - 10000 Hz with equation  B11
def B11(omega):
    """
    Equation B11 of [1].

    Parameters
    ----------
    omega : float
        angular frequency.

    Returns
    -------
    float
        damping loss of cavity.

    """
    # alpha = 1
    return air.c0*S_dwc_edge/6/omega/V_dwc

eta_DWC_data[0,14:] = B11(freq_tab_BI.data[14:])
# Create damping loss as Signal 
eta_DWC = mC.Signal(freq_tab_BI, eta_DWC_data, dof.DOF(0,0,dof.DOFtype(typeID = 1)))



# Crossmode Frequency f_d
f_d = air.c0/2/L2 
print('Crossmode frequency: {0:.1f}'.format(f_d))

#%% Double wall resonance
m1 = alu_3_18mm.mass_per_area
m2 = alu_6_36mm.mass_per_area
m3 = alu_0_88mm.mass_per_area
m4 = alu_9_26mm.mass_per_area

def DW_frequency (m1,m2,d):
    return 0.5/np.pi*np.sqrt(air.rho0*air.c0**2/d/m1/m2*(m1+m2))

print('Double wall resonance Fig.5: {0:.1f}Hz'.format(DW_frequency(m1,m1,L2)))
print('Double wall resonance Fig.6: {0:.1f}Hz'.format(DW_frequency(m1,m2,L2)))
print('Double wall resonance Fig.7: {0:.1f}Hz'.format(DW_frequency(m1,m3,L2)))



#%% Subsystem definition
# Plate subsystems
plate1_fig_5 = st2Dsys.RectangularPlate(2, L1,L3,prop=alu_3_18mm, eta = 0.005) 
plate2_fig_5 = st2Dsys.RectangularPlate(4, L1,L3,prop=alu_3_18mm, eta = 0.005)
 
plate1_fig_6 = st2Dsys.RectangularPlate(2, L1,L3,prop=alu_6_36mm, eta = 0.005) 
plate2_fig_6 = st2Dsys.RectangularPlate(4, L1,L3,prop=alu_3_18mm, eta = 0.005) 

plate1_fig_7 = st2Dsys.RectangularPlate(2, L1,L3,prop=alu_0_88mm, eta = 0.005) 
plate2_fig_7 = st2Dsys.RectangularPlate(4, L1,L3,prop=alu_0_88mm, eta = 0.005) 

plate1_fig_10 = st2Dsys.RectangularPlate(2, L1,L3,prop=alu_9_26mm, eta = 0.005) 
plate2_fig_10 = st2Dsys.RectangularPlate(4, L1,L3,prop=alu_9_26mm, eta = 0.005) 


# Cavity subsystems
room1 = ac3Dsys.Acoustic3DSystem(1, V1, 0, 0, air, absorption_area = As1, damping_type= ['surface'] )
room2 = ac3Dsys.Acoustic3DSystem(5, V2, 0, 0, air, absorption_area = As2, damping_type= ['surface'])

# Double wall cavity
# The flat_cavity_sw triggers the doubled radiation efficiency from plate to DWC
dwc     = ac3Dsys.RectangularRoom(3, L1, L2, L3, air, eta = eta_DWC,flat_cavity_sw = True)
# Extra cavity of Fig. 3 that uses other dimensions
dwc_fig3 = ac3Dsys.RectangularRoom(3, L1_fig3, L2, L3, air, eta = eta_DWC,flat_cavity_sw = True)

#%% Modal densities (of other modified cavity) Fig. 3
plt.figure(1)
plt.plot(freq,dwc_fig3.modal_density(om),label = 'Maa')
# Calculate modal density from mode counting
n_num,om_lim = dwc_fig3.modal_density_precise(om)

plt.plot(freq[1:],n_num, label = 'Mode counting')
#plt.plot(freq,dwc_fig3.modal_density(om,method = 'Price'),label = 'Price mod dens')
plt.plot(freq,(L1_fig3*L3*om/2/np.pi/air.c0**2),label = 'Eq (28)')
plt.plot(freq,(L1_fig3*L3*L2*om**2/2/np.pi**2/air.c0**3),label = 'Eq (29)')
#plt.plot(freq,dwc_fig3.modal_density(om,method = 'Price'),label = 'Price mod dens')

plt.annotate(r'$c_0/2L_2$', xy=(f_d,L1_fig3*L3*f_d/air.c0**2), xytext=(5000,0.05),
              arrowprops=aprop)
plt.title("Fig.3")
plt.xlim((100,1e4))
plt.ylim((1e-2,1))
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f/$Hz')
plt.ylabel('$n(\omega)/$s')
plt.xticks(fc,fclabels)
plt.legend()
plt.grid(which = 'both')
plt.show


#%% Junctions Fig 5 Model
J123_fig_5 = con.AreaJunction((room1,plate1_fig_5,dwc))
J345_fig_5 = con.AreaJunction((dwc,plate2_fig_5,room2))

#%% Junctions Fig 6 Model
J123_fig_6 = con.AreaJunction((room1,plate1_fig_6,dwc))
J345_fig_6 = con.AreaJunction((dwc,plate2_fig_6,room2))

#%% Junctions Fig 7 Model
J123_fig_7 = con.AreaJunction((room1,plate1_fig_7,dwc))
J345_fig_7 = con.AreaJunction((dwc,plate2_fig_7,room2))

#%% Junctions Fig 10 Model
J123_fig_10 = con.AreaJunction((room1,plate1_fig_10,dwc))
J345_fig_10 = con.AreaJunction((dwc,plate2_fig_10,room2))

#%% Load definition
# Source in cavity 1
power1mWatt = lC.Load(omega, 0.001*np.ones(omega.shape), 
                      dof.DOF(1,0,dof.DOFtype(typestr = 'power')), name = '1mWatt')


#%% Select to Plot Fig 5 (and 8),6 7 or 10
fig_nr = 10

if fig_nr == 5:
    # SEA Model used in Fig 5
    SEA_model = mds.HybridModel((room1,plate1_fig_5,dwc,plate2_fig_5,room2),xdata=omega)
    # Add junctions
    SEA_model.add_junction({'areaJ_123':J123_fig_5})
    SEA_model.add_junction({'areaJ_345':J345_fig_5})
elif fig_nr == 6:
    # SEA Model used in Fig 6
    SEA_model = mds.HybridModel((room1,plate1_fig_6,dwc,plate2_fig_5,room2),xdata=omega)
    # Add junctions
    SEA_model.add_junction({'areaJ_123':J123_fig_6})
    SEA_model.add_junction({'areaJ_345':J345_fig_6})
elif fig_nr == 7:
    # SEA Model used in Fig 6
    SEA_model = mds.HybridModel((room1,plate1_fig_7,dwc,plate2_fig_7,room2),xdata=omega)
    # Add junctions
    SEA_model.add_junction({'areaJ_123':J123_fig_7})
    SEA_model.add_junction({'areaJ_345':J345_fig_7})
elif fig_nr == 10:
    # SEA Model used in Fig 6
    SEA_model = mds.HybridModel((room1,plate1_fig_10,dwc,plate2_fig_10,room2),xdata=omega)
    # Add junctions
    SEA_model.add_junction({'areaJ_123':J123_fig_10})
    SEA_model.add_junction({'areaJ_345':J345_fig_10})
else:
    raise ValueError

# Add loads (same for all cases)
SEA_model.add_load('1mWatt',power1mWatt) # add 1Watt per band 
#%% Solve models Fig 5-7
if fig_nr in {5,6,7}:
    SEA_model.create_SEA_matrix(sym = 1)
    SEA_model.solve()

    # Derive paths of power input
    pow_in_room2 = SEA_model.power_input(5)

    #%% Engineering units velocity - plate results
    plt.figure(2)
    plt.title('Plate velocity of plates from fig {0}'.format(fig_nr))
    
    SEA_model.result.plot(2,ID=[2,4],dof=[3,3],xscale = 'log',yscale = 'log',
                           fulllegstr = ('plate B1','plate B2'),
                           xlabel = '$f/$Hz', ylabel = '$v_{\\rm rms}/$m s$^{-1}$',
                           xticks = 2*np.pi*fc[2:], xticklabels = fclabels[2:])
    
    
    #%% Engineering units pressure - cavity results
    plt.figure(3)
    plt.title('Cavity pressure from fig {0}'.format(fig_nr))
    SEA_model.result.plot(3,ID=[1,3,5],xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                        fulllegstr = ('room 1','DWC','room 2'),
                        xlabel = '$f/$Hz', ylabel = '$p_{\\rm rms}/$Pa',
                        xticks = 2*np.pi*fc[2:], xticklabels = fclabels[2:])
    
    
    
    #%% plot comparative modes in band 
    plt.figure(4)
    plt.plot(freq,plate1_fig_5.modes_in_band(om,3),label = 'plate 1 Fig 5')
    plt.plot(freq,plate1_fig_6.modes_in_band(om,3),label = 'plate 1 Fig 6')
    plt.plot(freq,plate1_fig_7.modes_in_band(om,3),label = 'plate 1 Fig 7')
    #plt.plot(freq,plate2.modes_in_band(om,3),label = 'plate 2')
    plt.plot(freq,dwc.modes_in_band(om),label = 'DWC')
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('$f/$Hz')
    plt.ylabel('N')
    plt.title('Modes in 3rd oct band')
    plt.ylim(bottom=0.1)
    plt.xticks(fc,fclabels)
    plt.legend()
    plt.grid(which = 'both')
    plt.show
 
#%% Plot Fig 9 radation efficiency
sigma1 = plate1_fig_5.radiation_efficiency(om,fluid = air)
sigma2 = plate1_fig_5.radiation_efficiency_simple(om,fluid = air)

# Show coincidence frequencies
f_c1 = alu_3_18mm.coincidence_frequency(air.c0)/2/np.pi
print('3.18 mm alu plate coincidence frequency = {0:.1f}'.format(f_c1))
f_c2 = alu_6_36mm.coincidence_frequency(air.c0)/2/np.pi
print('6.36 mm alu plate coincidence frequency = {0:.1f}'.format(f_c2))
f_c3 = alu_0_88mm.coincidence_frequency(air.c0)/2/np.pi
print('0.88 mm alu plate coincidence frequency = {0:.1f}'.format(f_c3))
f_c4 = alu_9_26mm.coincidence_frequency(air.c0)/2/np.pi
print('0.88 mm alu plate coincidence frequency = {0:.1f}'.format(f_c4))

plt.figure(5)
plt.plot(freq,10*np.log10(sigma1),'-.',label = 'Leppington')
plt.plot(freq,10*np.log10(sigma2),'-.',label = 'ISO EN 12354-1')
plt.title('Fig. 9')
plt.xlabel('$f/$Hz')
plt.ylabel('$\sigma/dB$')
plt.xscale('log')
plt.ylim(-50,25)
plt.xticks(fc,fclabels)
plt.legend()
plt.grid(which = 'both')


#%% Calculate transmission loss fig 5-7
# Calculate transmission loss from pressure difference
# Note: This is allowed because wall and absorption surface are equal S=As2
#%% Solve models Fig 5-7
if fig_nr in {5,6,7}:
    p1 = SEA_model.result[0].ydata.flatten()
    p2 = SEA_model.result[6].ydata.flatten()
    tau = (p2/p1)**2

#%% Plot transmission loss fig 5-7
    plt.figure(6)
    plt.semilogx(freq,-10*np.log10(tau),'-.',label = 'Pyva',color = 'C1')
    plt.xticks(fc,fclabels)
    plt.ylim((0,75))
    plt.xlabel('$f/$Hz')
    plt.ylabel('$TL/$dB')
    plt.title('Fig. {0}'.format(fig_nr))
    plt.grid(which = 'both')
    plt.legend()
    plt.show()

#%% Plot input power
    plt.figure(7)
    plt.title('Powert input into room 2 from fig {0}'.format(fig_nr))
    
    pow_in_room2.plot(7,yscale = 'log',fulllegstr = ['$\Pi_{DWC}$','$\Pi_{plate2}$'],
                      xscale = 'log',
                      xlabel = '$f/$Hz', ylabel = '$\Pi_{in}/$W',
                      xticks = 2*np.pi*fc[2:], xticklabels = fclabels[2:])

#%% Redo calculatiwith 5% damping if Fig 5 is selected to generate Fig 8
if fig_nr == 5: # Initiates Fig 8 calulation 
    plate1_fig_5.eta = 0.05
    plate2_fig_5.eta = 0.05

    SEA_model.create_SEA_matrix(sym = 1)
    SEA_model.solve()  

    p1 = SEA_model.result[0].ydata.flatten()
    p2 = SEA_model.result[6].ydata.flatten()
    tau_eta_0_05 = (p2/p1)**2
 
    plt.figure(8)
    plt.semilogx(freq,-10*np.log10(tau),'-.',label = r'$\eta = 0.5$%', color = 'C1')
    plt.semilogx(freq,-10*np.log10(tau_eta_0_05),'-.',label = r'$\eta = 5$%', color = 'C2')
    plt.xticks(fc[2:],fclabels[2:])
    plt.ylim((0,75))
    plt.xlabel('$f/$Hz')
    plt.ylabel('$TL/$dB')
    plt.title('Fig 8')
    plt.grid(which = 'both')
    plt.legend()
    plt.show()


#%% calculate DWC thickness variation of Fig 10 

thicknesses = [0.01,0.2,0.4]
tau = [None]*3
if fig_nr == 10:
    plt.figure(6)

    for i_cav,t in enumerate(thicknesses):
        print('Evaluating thickness {0:.1f}mm'.format(t*1000))
        dwc.Ly = t
        dwc.update()
        print('Double wall frequency for L2={0:.0f}cm : {1:.1f}Hz'.format(t*100,DW_frequency(m4,m4,t)))
        
        SEA_model.create_SEA_matrix(sym = 1)
        SEA_model.solve()  

        p1 = SEA_model.result[0].ydata.flatten()
        p2 = SEA_model.result[6].ydata.flatten()
        tau[i_cav] = (p2/p1)**2
        
        plt.semilogx(freq,-10*np.log10(tau[i_cav])+0*i_cav,label = '$L_2=${0:.0f}cm'.format(t*100))
    
    plt.xticks(fc,fclabels)
    plt.ylim((0,75))
    plt.xlabel('$f/$Hz')
    plt.ylabel('$TL/$dB')
    plt.title('Fig. {0}'.format(fig_nr))
    plt.grid(which = 'both')
    plt.legend()
    plt.show()
    

#%% Plot damping
#%% Engineering units pressure - cavity results
plt.figure(7)
plt.title('Double wall cavity damping')
plt.grid(which='both')

eta_DWC.plot(7,xscale = 'log', fulllegstr = ('$\eta_3$',),
                    xlabel = '$f/$Hz', ylabel = r'$\eta$',
                    xticks = 2*np.pi*fc, xticklabels = fclabels)




#%% Create DW cavity and plate2 as trim
air_damped = matC.Fluid(eta = eta_DWC)
il_air_7_11cm = iL.FluidLayer(L2,air_damped)
il_plate2     = iL.PlateLayer(alu_3_18mm)