% Calibrated simulation with observed data — soil and routing parameters (set manually).

clc;
clear;
close all;
addpath("modules/parameters", "modules/model");

%% PARAMETERS

shared = struct( ...
    'baseflow_mm_h', 0.0, ...   % [mm/h] 
    's0', 0.484, ...
    'LQ_k', NaN, ...            % NaN = default Valfredda
    'LQ_theta', NaN);

gamCase = struct( ...
    'label', 'Gamma', ...
    'Kc', 0.83, ...             % [-]
    'b', 1.8, ...               % [-]
    'Ks', 101.7, ...            % [mm/h]
    'k', 0.018, ...             % [1/h]
    'uh', 20.2, ...             % [m/h]
    'beta', 160, ...            % [-] par.alpha.beta
    'theta', 0.745, ...         % forma IUH a Lmax [-]
    'LQ_k', NaN, ...
    'LQ_theta', NaN);

%% Setup 
par0 = makeDefaultParameters("observed", "Valfredda");
par0.domainT.dt = 1;
bf = shared.baseflow_mm_h;

fprintf('=== calibratedSimulation ===\n\n');
printCaseParams(gamCase);

%% Simulations (soil + routing)
[rt, par, ~] = runGammaCase(par0, gamCase, shared);

[Q_obs, obsMeta] = readObservedQ(par0);

NtSim = height(rt);
NtObs = numel(Q_obs);
stepsPerDay = round(24 / par0.domainT.dt);

%% Zoom 15--19
if ismember('date', rt.Properties.VariableNames)
    datesSim = datetime(rt.date, 'InputFormat', 'preserve');
else
    meteoH = readHourlyMeteo(par0, false);
    datesSim = datetime(meteoH.date(1:min(height(meteoH), NtSim)), 'InputFormat', 'preserve');
end
iMay = find(isMay(datesSim));
nCmp = min(numel(iMay), NtObs);
jZoom = find(day(datesSim(iMay)) >= 15 & day(datesSim(iMay)) <= 19);
iZoom = iMay(jZoom);
iZoom = iZoom(iZoom >= 1 & iZoom <= NtSim);
jZoomCmp = jZoom(jZoom <= nCmp);
t_zoom = day(datesSim(iZoom)) + hour(datesSim(iZoom)) / 24 ...
    + minute(datesSim(iZoom)) / 1440;

%% Error Metrics 15--19 maggio
hEvt = hydroMetrics(rt.Q(iZoom(1:numel(jZoomCmp))), Q_obs(jZoomCmp));
fprintf('\n--- Error metrics ---\n');
fprintf('NSE=%.3f | KGE=%.3f | RMSE=%.3f | score=%.3f\n', ...
    hEvt.NSE, hEvt.KGE, hEvt.RMSE, 0.5*hEvt.NSE + 0.5*hEvt.KGE);

%% Colors 
cObs = [0.65, 0.85, 0.95];
cSim = [0.10, 0.20, 0.40];
cL = [0.25, 0.45, 0.75];      
cP = [0.45 0.55 0.70];
cQs = [0.25, 0.60, 0.40];
cQd = [0.35, 0.55, 0.75];
cQtot = [0.85, 0.65, 0.70];
% Dimensions 
figW = 9;    % cm
figH = 16;   % cm
figPos = [3, 2, figW, figH];

%% PLOT (zoom 15--19 May)
figure('Name', 'Gamma — P | L | Q | componenti (15--19 mag)', 'NumberTitle', 'off', 'Color', 'w', ...
    'WindowStyle', 'normal', 'Units', 'centimeters', 'Position', figPos, ...
    'PaperUnits', 'centimeters', 'PaperSize', [figW figH], ...
    'PaperPosition', [0 0 figW figH]);

% characters dimensions
tl = tiledlayout(4, 1, 'TileSpacing', 'compact', 'Padding', 'compact');
fsLabel = 12;   % axes
fsLegend = 10;  % legend
fsTick = 9;     % tick numbers
fsPanel = 10;   % letters
xTickMay = 15:19;
xTickLbl = {'$\mathrm{May\,15}$', '$\mathrm{May\,16}$', '$\mathrm{May\,17}$', ...
    '$\mathrm{May\,18}$', '$\mathrm{May\,19}$'};

% Rain
ax1 = nexttile(tl);
bar(t_zoom, rt.P(iZoom), 1, 'FaceColor', cP, 'EdgeColor', 'none');
set(ax1, 'YDir', 'reverse', 'XTickLabel', [], 'TickLabelInterpreter', 'latex', ...
    'FontSize', fsTick, 'Box', 'off', 'TickDir', 'in', 'TickLength', [0.015 0.025]);
ylabel('$h \,[mm/h]$', 'Interpreter', 'latex', 'FontSize', fsLabel);
xlim([15 19]);
xticks(xTickMay);
addPanelLetter(ax1, 'a)', fsPanel);
drawTopRightBox(ax1);

%  L(t) [km] 
ax2 = nexttile(tl);
hold on;
plot(t_zoom, rt.L(iZoom) / 1000, '-', 'Color', cL, 'LineWidth', 1.4);
set(ax2, 'XTickLabel', [], 'TickLabelInterpreter', 'latex', 'FontSize', fsTick, ...
    'Box', 'off', 'TickDir', 'in', 'TickLength', [0.015 0.025]);
ylabel('$L\, [km]$', 'Interpreter', 'latex', 'FontSize', fsLabel);
xlim([15 19]);
xticks(xTickMay);
addPanelLetter(ax2, 'b)', fsPanel);
drawTopRightBox(ax2);

% Discharge
ax3 = nexttile(tl);
hold on;
plot(t_zoom, rt.Q(iZoom), '-', 'Color', cSim, 'LineWidth', 1.4, 'DisplayName', '$dynamic\, \textit{IUH}\, model$');
plot(t_zoom(1:numel(jZoomCmp)), Q_obs(jZoomCmp), '.', 'Color', cObs, 'MarkerSize', 10, ...
    'DisplayName', '$\textit{DATA}$');
set(ax3, 'XTickLabel', [], 'TickLabelInterpreter', 'latex', 'FontSize', fsTick, ...
    'Box', 'off', 'TickDir', 'in', 'TickLength', [0.015 0.025]);
ylabel('$Q\, [mm/h]$', 'Interpreter', 'latex', 'FontSize', fsLabel);
xlim([15 19]);
xticks(xTickMay);
ylim([0 2]);
yticks(0:0.5:2);
lgd3 = legend(ax3, 'Location', 'northwest', 'Box', 'off', 'Interpreter', 'latex', 'FontSize', 8);
lgd3.ItemTokenSize = [15, 9];   % linea blu dimezzata rispetto al default (~30)
addPanelLetter(ax3, 'c)', fsPanel);
drawTopRightBox(ax3);

% Qs (stat), Qd (dyn) 
ax4 = nexttile(tl);
hold on;
plot(t_zoom, rt.Qs(iZoom), ':', 'Color', cQs, 'LineWidth', 1.3, 'DisplayName', '$Q_{st}$');
plot(t_zoom, rt.Qd(iZoom), '--', 'Color', cQd, 'LineWidth', 1.2, 'DisplayName', '$Q_{dyn}$');
set(ax4, 'TickLabelInterpreter', 'latex', 'FontSize', fsTick, ...
    'Box', 'off', 'TickDir', 'in', 'TickLength', [0.015 0.025]);
xlabel('$2024$', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('$Q_{st}, Q_{dyn} \, [mm/h]$', 'Interpreter', 'latex', 'FontSize', fsLabel);
xlim([15 19]);
xticks(xTickMay);
xticklabels(xTickLbl);
lgd4 = legend(ax4, 'Location', 'east', 'Box', 'off', 'Interpreter', 'latex', 'FontSize', 10);
styleLegend(lgd4, fsLegend);
% nord-est ma a meta' dell'asse y: sposta un filo verso l'alto dal centro destro
ax4.Units = 'normalized';
lgd4.Units = 'normalized';
pAx = ax4.Position;
pLg = lgd4.Position;
lgd4.Position = [pAx(1) + pAx(3) - pLg(3) - 0.01, ...
                 pAx(2) + 0.55*pAx(4) - 0.5*pLg(4), ...
                 pLg(3), pLg(4)];
addPanelLetter(ax4, 'd)', fsPanel);
drawTopRightBox(ax4);

linkaxes([ax1, ax2, ax3, ax4], 'x');

%% FUNCTIONS

function [rt, par, soil] = runGammaCase(par0, casePar, shared)
    par = par0;
    par.rain.soil.s0 = shared.s0;
    par.rain.et.Kc = casePar.Kc;
    par.rain.soil.b = casePar.b;
    par.rain.soil.Ks = casePar.Ks;
    par.f0.k = casePar.k;
    par.f0.uh = casePar.uh;
    par.alpha.beta = casePar.beta;
    par.f0.theta = casePar.theta;
    lqK = pickParam(casePar, shared, 'LQ_k');
    lqTh = pickParam(casePar, shared, 'LQ_theta');
    if ~isnan(lqK), par.LQ.k = lqK; end
    if ~isnan(lqTh), par.LQ.theta = lqTh; end
    par = makeDependentParameters(par, 'rain.baseflow');
    par.rain.baseflow = shared.baseflow_mm_h;

    soil = par.rain.type(par.domainT, par.rain);
    F_L_tau = analyticalSolution(par.domainL.Lmin, par.domainL.Lmax, par.domainL.dL, ...
        par.domainTau.tau, par.f0.uh, par.alpha.beta, par.f0.k, par.f0.theta);
    rt = soil;
    [rt, ~] = simulateDynamic(rt, par.alpha.fun, par.f0.uh, par.LQ, F_L_tau, par.rain.baseflow);
end

function v = pickParam(casePar, shared, name)
    if isfield(casePar, name) && ~isnan(casePar.(name))
        v = casePar.(name);
    elseif isfield(shared, name) && ~isnan(shared.(name))
        v = shared.(name);
    else
        v = NaN;
    end
end

function printCaseParams(casePar)
    fprintf('  Soil:  Kc=%.3f [-] | b=%.3f [-] | Ks=%.1f mm/h\n', casePar.Kc, casePar.b, casePar.Ks);
    fprintf('  Routing: k=%.5f [1/h] | uh=%.1f [m/h]| beta=%.1f [h/m] | f0.theta=%.3f\n', ...
        casePar.k, casePar.uh, casePar.beta, casePar.theta);
    if isfield(casePar, 'LQ_k') && ~isnan(casePar.LQ_k)
        fprintf('  L(Q): k=%.3f | theta=%.4f mm/h\n', casePar.LQ_k, casePar.LQ_theta);
    end
end

function tf = isMay(d)
    d = datetime(d, 'InputFormat', 'preserve');
    tf = month(d) == 5 & day(d) >= 1 & day(d) <= 31;
end

function meteoH = readHourlyMeteo(par, useFull)
    if nargin < 2, useFull = false; end
    meteoRaw = readtable(par.rain.fileInput);
    spdModel = round(24 / par.domainT.dt);
    meteoH = aggregateToModelDt(meteoRaw, par.rain.spdFile, spdModel, 'meteo');
    if ~useFull
        Nt = min(par.domainT.Nt, height(meteoH));
        meteoH = meteoH(1:Nt, :);
    end
end

function drawTopRightBox(ax)
    xl = ax.XLim;
    yl = ax.YLim;
    if strcmp(ax.YDir, 'reverse')
        yTop = yl(1);
    else
        yTop = yl(2);
    end
    hold(ax, 'on');
    plot(ax, [xl(1) xl(2)], [yTop yTop], 'k-', 'LineWidth', ax.LineWidth, ...
        'HandleVisibility', 'off', 'Clipping', 'off');
    plot(ax, [xl(2) xl(2)], yl, 'k-', 'LineWidth', ax.LineWidth, ...
        'HandleVisibility', 'off', 'Clipping', 'off');
end

function styleLegend(lgd, fs)
    lgd.FontSize = fs;
    lgd.Interpreter = 'latex';
    lgd.Box = 'off';
    if isprop(lgd, 'ItemTokenSize')
        lgd.ItemTokenSize = [20, 9];
    end
end

function addPanelLetter(ax, letter, fs)
% Font "normale" sans-serif bold (come nello screen), non LaTeX
    text(ax, 0.97, 0.97, letter, ...
        'Units', 'normalized', ...
        'HorizontalAlignment', 'right', 'VerticalAlignment', 'top', ...
        'Interpreter', 'none', 'FontName', 'Helvetica', ...
        'FontWeight', 'bold', 'FontSize', fs);
end
