clear
warning('off')
format compact
load Case1_Chl1_HGg0.95_theta_0.mat

N = 100000; % photons
gamma_cutoff = 10; % angle for transport approximation
theta_s = 0; % sun position
theta_0 = asind(sind(theta_s)/1.34); % incident angle. After refraction theta_0 = asind(sind(theta_s)/1.34); 60->40

%%
ii = 1;
N_reduction = 4; % calc every N-th wl
for i = 1:N_reduction:length(a)
    [Kd_layers(ii), R0(ii)] = MC_3D_2024_old(a(i),b(i),bw(i),theta_s,N,gamma_cutoff);
    [num2str(100*ii*N_reduction/length(a)) '%']
    ii = ii+1;
end

%% Gordon Kd
% calculate bb from joint VSF
P_water = @(gamma) 0.06083*(1+0.92493*(cosd(gamma)).^2); % water

g = 0.95; % HG Scattering Phase Function
gamma_values = [0:0.001:180]';
P = @(gamma) (1-g.^2)./(1+g.^2-2*g*cosd(gamma)).^(3/2)/4/pi;

beta_int = (P_water(gamma_values)*bw' + P(gamma_values)*(b-bw)')./b'; % VSF_integral = 1
cos_th = cosd(gamma_values);
d_cos = diff([1; cos_th]);

% VSF_integral = -2*pi*sum(beta_int.*d_cos) % should be equal to 1. '-' from inversion of limits 
% g_test = -2*pi*sum(beta_int.*d_cos.*cos_th) % % should be equal to g

ind = gamma_values>90;
bb_b_ratio = -2*pi*(sum(beta_int(ind,:).*d_cos(ind)))';
bb = b.*bb_b_ratio;

% Gordon'89 Kd
g_diff_fraction = 0; % 0 for the Sun in the black sky
D0 = (1-g_diff_fraction)/cosd(theta_0)+1.197*g_diff_fraction;

F = 1 - bb./b;
w = b./(a+b);

Kd_Gordon_deep = (a+b).*D0.*(1.3197*(1-w.*F) - 0.7559*(1-w.*F).^2 + 0.4655*(1-w.*F).^3); % in z(Ed=0.1Ed(0)) layer
% Kd_Gordon_surf = 1.04*D0*(a + bb); % under the surface

clearvars P P_water gamma_values beta_int cos_th D0 F d_cos g_diff_fraction

%% true HL Kd
Kd_HL_true = zeros(1, size(Ed_hl,1)-1);
for i = 1:size(Ed_hl,1)-1
    Ed_hl_calc = Ed_hl(i+1,Ed_hl(i+1,:)>Ed_hl(i+1,1)/10); % in z_0.1 layer, like in Gordon
    Ed_fit = robustfit(Ed_hl(1,Ed_hl(i+1,:)>Ed_hl(i+1,1)/10), log(Ed_hl_calc)); % bisquare
    Kd_HL_true(i) = abs(Ed_fit(2));
end
clearvars Ed_hl_calc Ed_fit

%% plots
clf
wl = 400:20*N_reduction:700;
wl_Hl = 400:20:700;

subplot(1,2,1)
% plot(wl_Hl, Kd_HL_surf, 'r', 'LineWidth', 2)
hold on, grid on
plot(wl_Hl, Kd_HL_true, 'k', 'LineWidth', 3)
% plot(wl_Hl, Kd_Gordon_surf, 'g', 'LineWidth', 2)
plot(wl_Hl, Kd_Gordon_deep, 'g--', 'LineWidth', 3)
plot(wl, Kd_layers, 'o', 'LineWidth', 3)

legend('HydroLight', 'Gordon', 'MC', 'Location', 'northwest')
xlabel('\lambda, íì')
ylabel('K_d, ì^{-1}')
set(gca, 'LineWidth', 1.5, 'FontSize', 14)

subplot(1,2,2)
plot(wl_Hl, R, 'k', 'LineWidth', 3)
hold on, grid on
plot(wl, R0, 'o', 'LineWidth', 3)
legend('HydroLight', 'MC')

xlabel('\lambda, íì')
ylabel('R')
set(gca, 'LineWidth', 1.5, 'FontSize', 14)
