%%
clear; close all; clc;

%%
load('./mat/Gd_ol_cl.mat', 'Gd_ol_20', 'Gd_cl_20');

%%
load('./mat/perturbations.mat', 'Wxg')

%%
load('./mat/G_gm_to_dh.mat', 'G_gm_to_dh')
load('./mat/psd_ground_motion.mat', 'psd_f', 'pxx')

%%
bodeFig({Gd_ol_20(1, 1), G_gm_to_dh})

%% PSD
freqs = logspace(-2, 2, 1000);

gm_ol = abs(squeeze(freqresp(Wxg*Gd_ol_20(1, 1), freqs, 'Hz')));
gm_cl = abs(squeeze(freqresp(Wxg*Gd_cl_20(1, 1), freqs, 'Hz')));

figure;
hold on;
plot(freqs, gm_ol)
plot(freqs, gm_cl)
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [$Hz$]'); ylabel('PSD [$m/\sqrt{Hz}$]');
hold off;

%% CAS
freqs = logspace(-1, 2, 1000);

gm_ol = abs(squeeze(freqresp(Wxg*Gd_ol_20(1, 1), freqs, 'Hz')));
gm_cl = abs(squeeze(freqresp(Wxg*Gd_cl_20(1, 1), freqs, 'Hz')));

dw = freqs - [0, freqs(1:end-1)];

figure;
hold on;
plot(freqs, cumsum(gm_ol'.*dw))
plot(freqs, cumsum(gm_cl'.*dw))
set(gca, 'XScale', 'log');
% set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('CAS [m]');
hold off;