<reponame>antonyvm1102/Redbox<filename>plot_data_RMS.py<gh_stars>0
import matplotlib.pyplot as plt
import numpy as np

"""
Plot RMS for x-/y-/z-signal vs 1/3 octave band frequencies and compare to VC curves.
"""

# rms_x_all = np.loadtxt("14208_betacampus_pos1_rms_x_all.txt")
# rms_y_all = np.loadtxt("14208_betacampus_pos1_rms_y_all.txt")
# rms_z_all = np.loadtxt("14208_betacampus_pos1_rms_z_all.txt")
# rms_x_all = np.loadtxt("14208_betacampus_pos2_rms_x_all.txt")
# rms_y_all = np.loadtxt("14208_betacampus_pos2_rms_y_all.txt")
# rms_z_all = np.loadtxt("14208_betacampus_pos2_rms_z_all.txt")
# rms_x_all = np.loadtxt("14208_Huygensgebouw_rms_x_all.txt")
# rms_y_all = np.loadtxt("14208_Huygensgebouw_rms_y_all.txt")
# rms_z_all = np.loadtxt("14208_Huygensgebouw_rms_z_all.txt")
# rms_x_all = np.loadtxt("14208_betacampus_pos1_rms_x_weekday.txt")
# rms_y_all = np.loadtxt("14208_betacampus_pos1_rms_y_weekday.txt")
# rms_z_all = np.loadtxt("14208_betacampus_pos1_rms_z_weekday.txt")
# rms_x_all = np.loadtxt("14208_betacampus_pos2_rms_x_weekday.txt")
# rms_y_all = np.loadtxt("14208_betacampus_pos2_rms_y_weekday.txt")
# rms_z_all = np.loadtxt("14208_betacampus_pos2_rms_z_weekday.txt")
rms_x_all = np.loadtxt("14208_Huygensgebouw_rms_x_weekday.txt")
rms_y_all = np.loadtxt("14208_Huygensgebouw_rms_y_weekday.txt")
rms_z_all = np.loadtxt("14208_Huygensgebouw_rms_z_weekday.txt")

f_band = rms_x_all[..., 0]

def VC_curve(v_max, f_trans, f_band):
    v = np.zeros_like(f_band)
    for i in range(len(f_band)):
        if f_band[i] < f_trans:
            v[i] = v_max * f_trans / f_band[i]
        else:
            v[i] = v_max
    return v

v_VCA = VC_curve(50./1000., 8., f_band)
v_VCB = VC_curve(25./1000., 8., f_band)
v_VCC = VC_curve(12.5/1000., 8., f_band)
v_VCD = VC_curve(6.25/1000., 8., f_band)
v_VCE = VC_curve(3.125/1000., 8., f_band)
v_VCF = VC_curve(1.5625/1000., 8., f_band)
v_VCG = VC_curve(0.78125/1000., 8., f_band)
v_VCH = VC_curve(0.390625/1000., 8., f_band)

plt.subplot(3,1,1)
plt.plot(rms_x_all[..., 0], rms_x_all[..., 5], c = 'k', ls = 'dotted', label = 'maximum')
# plt.plot(rms_x_all[..., 0], rms_x_all[..., 4], c = 'k', ls = 'dashed', label = 'mean + st.dev.')
plt.plot(rms_x_all[..., 0], rms_x_all[..., 3], c = 'k', ls = 'solid', label = 'mean')
# plt.plot(rms_x_all[..., 0], rms_x_all[..., 2], c = 'k', ls = 'dashed', label = 'mean - st.dev.')
plt.plot(rms_x_all[..., 0], rms_x_all[..., 1], c = 'k', ls = 'dotted', label = 'minimum')
plt.fill_between(rms_x_all[..., 0], rms_x_all[... , 1], rms_x_all[..., 5], color = 'lightgray')
plt.plot(f_band, v_VCA, '-', label = 'VC-A')
plt.plot(f_band, v_VCB, '-', label = 'VC-B')
plt.plot(f_band, v_VCC, '-', label = 'VC-C')
plt.plot(f_band, v_VCD, '-', label = 'VC-D')
plt.plot(f_band, v_VCE, '-', label = 'VC-E')
plt.plot(f_band, v_VCF, '-', label = 'VC-F')
plt.plot(f_band, v_VCG, '-', label = 'VC-G')
plt.plot(f_band, v_VCH, '-', label = 'VC-H')
plt.legend(bbox_to_anchor=(-0.05,-0.5))
plt.xscale('log')
plt.xlim([np.power(10,-0.2), 100])
plt.yscale('log')
plt.grid(which = 'both', lw = 0.3)
plt.xlabel('1/3 octave band frequency [Hz]')
plt.ylabel('RMS of x-signal [mm/s]')
plt.subplot(3,1,2)
plt.plot(rms_y_all[..., 0], rms_y_all[..., 5], c = 'k', ls = 'dotted', label = 'maximum')
# plt.plot(rms_y_all[..., 0], rms_y_all[..., 4], c = 'k', ls = 'dashed', label = 'mean + st.dev.')
plt.plot(rms_y_all[..., 0], rms_y_all[..., 3], c = 'k', ls = 'solid', label = 'mean')
# plt.plot(rms_y_all[..., 0], rms_y_all[..., 2], c = 'k', ls = 'dashed', label = 'mean - st.dev.')
plt.plot(rms_y_all[..., 0], rms_y_all[..., 1], c = 'k', ls = 'dotted', label = 'minimum')
plt.fill_between(rms_y_all[..., 0], rms_y_all[... , 1], rms_y_all[..., 5], color = 'lightgray')
plt.plot(f_band, v_VCA, '-', label = 'VC-A')
plt.plot(f_band, v_VCB, '-', label = 'VC-B')
plt.plot(f_band, v_VCC, '-', label = 'VC-C')
plt.plot(f_band, v_VCD, '-', label = 'VC-D')
plt.plot(f_band, v_VCE, '-', label = 'VC-E')
plt.plot(f_band, v_VCF, '-', label = 'VC-F')
plt.plot(f_band, v_VCG, '-', label = 'VC-G')
plt.plot(f_band, v_VCH, '-', label = 'VC-H')
plt.xscale('log')
plt.xlim([np.power(10,-0.2), 100])
plt.yscale('log')
plt.grid(which = 'both', lw = 0.3)
plt.xlabel('1/3 octave band frequency [Hz]')
plt.ylabel('RMS of y-signal [mm/s]')
plt.subplot(3,1,3)
plt.plot(rms_z_all[..., 0], rms_z_all[..., 5], c = 'k', ls = 'dotted', label = 'maximum')
# plt.plot(rms_z_all[..., 0], rms_z_all[..., 4], c = 'k', ls = 'dashed', label = 'mean + st.dev.')
plt.plot(rms_z_all[..., 0], rms_z_all[..., 3], c = 'k', ls = 'solid', label = 'mean')
# plt.plot(rms_z_all[..., 0], rms_z_all[..., 2], c = 'k', ls = 'dashed', label = 'mean - st.dev.')
plt.plot(rms_z_all[..., 0], rms_z_all[..., 1], c = 'k', ls = 'dotted', label = 'minimum')
plt.fill_between(rms_z_all[..., 0], rms_z_all[... , 1], rms_z_all[..., 5], color = 'lightgray')
plt.plot(f_band, v_VCA, '-', label = 'VC-A')
plt.plot(f_band, v_VCB, '-', label = 'VC-B')
plt.plot(f_band, v_VCC, '-', label = 'VC-C')
plt.plot(f_band, v_VCD, '-', label = 'VC-D')
plt.plot(f_band, v_VCE, '-', label = 'VC-E')
plt.plot(f_band, v_VCF, '-', label = 'VC-F')
plt.plot(f_band, v_VCG, '-', label = 'VC-G')
plt.plot(f_band, v_VCH, '-', label = 'VC-H')
plt.xscale('log')
plt.xlim([np.power(10,-0.2), 100])
plt.yscale('log')
plt.grid(which = 'both', lw = 0.3)
plt.xlabel('1/3 octave band frequency [Hz]')
plt.ylabel('RMS of z-signal [mm/s]')
plt.plot()
plt.subplots_adjust(top=0.95, bottom=0.05, left=0.15, right=0.95, hspace=0.3)
# plt.suptitle('Effective velocities at Betacampus position 1 during weekdays between 6:00 and 18:00 hours')
# plt.suptitle('Effective velocities at Betacampus position 2 during weekdays between 6:00 and 18:00 hours')
plt.suptitle('Effective velocities at Huygensgebouw during weekdays between 6:00 and 18:00 hours')
plt.show()