%% Calculate Rectangular Bay resonance frequency using Sutherland & Garrett (2005) analytical model

clear

load complex_A_infer.mat
AAh = AA(17,:);
AAm = AA(5,:);
load complex_A_10m_infer.mat
AAh_mhk = AA(17,:);
AAm_mhk = AA(5,:);

nh = {'M2','S2','N2','K2','K1','O1','P1','Q1'};     %harmonics to apply

% Numerical model output
Ag = AAh./AAm;                        %all constituent complex elevation gain
Ag_mhk = AAh_mhk./AAm_mhk;                        

% all constituent non-complex elevation gain
Arg = abs(Ag);
Arg_r = real(Ag);
Arg_i = imag(Ag);
Arg_mhk = abs(Ag_mhk);
Arg_r_mhk = real(Ag_mhk);
Arg_i_mhk = imag(Ag_mhk);

% all constituent non-complex phase normalized
pd = angle(Ag);
pd_mhk = angle(Ag_mhk);

%Analytical model. Vary resonant frequency, omega_0, and quadratic
%friction, lambda, to find best fit to model output
lambda = linspace(1e-6,1e-4,500);
omega_0 = linspace(5.82e-5,3.49e-4,500);    %varying from 10 to 30 h

for i = 1:length(lambda)
    for j = 1:length(omega_0)        
        for k = 1:length(nh)
                phi = (omega(k)/omega_0(j))*(pi/2);
                ep = (lambda(i)/(2*omega_0(j)))*(pi/2);
                kL = phi + 1i*ep;
                Ag_ana(i,j,k) = sec(kL);
            
            %Find minimum difference btw model output and analytical model
            diff(i,j,k) = abs(Ag(k)-Ag_ana(i,j,k)).^2;
            diff_mhk(i,j,k) = abs(Ag_mhk(k)-Ag_ana(i,j,k)).^2;
        end        
    end
end

% Sum the differences over all harmonics and locate minimum then
% corresponding lambda, omega_0 values
diff_s = sum(diff,3,'omitnan');
diff_s_mhk = sum(diff_mhk,3,'omitnan');

min_diff = min(diff_s,[],'all');
min_diff_mhk = min(diff_s_mhk,[],'all');

[i_lambda,i_omega_0] = find(diff_s == min_diff);
[i_lambda_mhk,i_omega_0_mhk] = find(diff_s_mhk == min_diff_mhk);

lambda_f = lambda(i_lambda);
lambda_f_mhk = lambda(i_lambda_mhk);

omega_0_f = omega_0(i_omega_0);
omega_0_f_mhk = omega_0(i_omega_0_mhk);

T0_f = (1/(omega_0_f/(2*pi)))/3600;
T0_f_mhk = (1/(omega_0_f_mhk/(2*pi)))/3600;

% Calculate quality factor, Q, for best fit
Q = omega_0_f/lambda_f;
Q_M2 = 1.5*Q;
Q_mhk = omega_0_f_mhk/lambda_f_mhk;
Q_M2_mhk = 1.5*Q_mhk;

% Calculate analytical model result for range of omega (can use same range
% used in fitting) and the best-fit lambda
om_a = omega_0;

phi = (om_a./omega_0_f).*(pi/2);
ep = (lambda_f/(2*omega_0_f))*(pi/2);
kL = phi + 1i*ep;
Ag_ana_f = sec(kL);

phi_mhk = (om_a./omega_0_f_mhk).*(pi/2);
ep_mhk = (lambda_f_mhk/(2*omega_0_f_mhk))*(pi/2);
kL_mhk = phi_mhk + 1i*ep_mhk;
Ag_ana_f_mhk = sec(kL_mhk);

Ag_ana_a_f = abs(Ag_ana_f);
Ag_ana_a_f_mhk = abs(Ag_ana_f_mhk);

Ag_ana_r_f = real(Ag_ana_f);
Ag_ana_r_f_mhk = real(Ag_ana_f_mhk);

Ag_ana_i_f = imag(Ag_ana_f);
Ag_ana_i_f_mhk = imag(Ag_ana_f_mhk);

ph_ana_n = angle(Ag_ana_f);
ph_ana_n_mhk = angle(Ag_ana_f_mhk);

% Sensitivity of constituents to friction / resonance
om_M2 = abs(1-omega(1)/omega_0_f);
om_M2_mhk = abs(1-omega(1)/omega_0_f_mhk);

om_K1 = abs(1-omega(5)/omega_0_f);
om_K1_mhk = abs(1-omega(5)/omega_0_f_mhk);

Qi_M2 = 0.5*(1/Q_M2);
Qi_M2_mhk = 0.5*(1/Q_M2_mhk);

Qi_K1 = 0.5*(1/Q);
Qi_K1_mhk = 0.5*(1/Q_mhk);


%% Compare mhk to no mhk cases

figure
subplot(211)
grid on
hold on
plot((1./(om_a./(2.*pi)))/3600,Ag_ana_a_f,'k','linewidth',2)
plot((1./(om_a./(2.*pi)))/3600,Ag_ana_a_f_mhk,'r','linewidth',2)
scatter((1./(omega./(2.*pi)))/3600,Arg,'k','linewidth',1.5)
scatter((1./(omega./(2.*pi)))/3600,Arg_mhk,'r','linewidth',1.5)
xlim([8 28])
% ylabel('Amplitude Gain Ratio')
% title('Admiralty Inlet - South Sound')
% legend('Base Case','Turbine Case')
set(gca,'FontSize',14,'XTickLabels',[])

subplot(212)
grid on
hold on
plot((1./(om_a./(2.*pi)))/3600,rad2deg(ph_ana_n),'k','linewidth',2)
plot((1./(om_a./(2.*pi)))/3600,rad2deg(ph_ana_n_mhk),'r','linewidth',2)
scatter((1./(omega./(2.*pi)))/3600,rad2deg(pd),'k','linewidth',1.5)
scatter((1./(omega./(2.*pi)))/3600,rad2deg(pd_mhk),'r','linewidth',1.5)
ylim([0 150])
xlim([8 28])
% ylabel('Phase Difference [deg.]')
% xlabel('Period, T [h.]')
set(gca,'FontSize',14)