function result_t = generateObservedRainfall(domainT, rain)
% GENERATEOBSERVEDRAINFALL 
%
% Legge rain.fileInput (tabella meteo sub-giornaliera o giornaliera),
% aggrega a domainT.dt [h], stima ETmax/ETsw da meteo, chiama soilMoistureDynamics.
%
% Output (table): t [h], P [mm], h [mm] leaching (P_sub), s [-], L [mm], EVT [mm].


    meteo = readMeteoFile(rain.fileInput);
    spdFile = getRainParam(rain, 'spdFile', 144);
    spdModel = round(24 / domainT.dt);
    meteo = aggregateToModelDt(meteo, spdFile, spdModel, 'meteo');

    P_mm = meteo.PRCP(:);
    nAvail = numel(P_mm);
    Nt = min(domainT.Nt, nAvail);
    % if nAvail < domainT.Nt
    %     warning('generateObservedRainfall: meteo (%d passi) piu'' corto di domainT.Nt (%d).', nAvail, domainT.Nt);
    % end

    P_mm = P_mm(1:Nt);
    meteo = meteo(1:Nt, :);
    t = domainT.t(1:Nt)';

    dt_h = domainT.dt;
    [ETmax, ETsw] = estimateSoilET(meteo, rain, dt_h);

    soil = prepareSoilForMmH(rain.soil);

    ensurePlotColors();
    [s, L, EVT] = soilMoistureDynamics(P_mm, dt_h, soil.s0, soil.sw, soil.sstar, ...
        soil.Ks, soil.b, ETsw, ETmax, soil.n, soil.Zr);

    result_t = table();
    result_t.t = t;
    result_t.P = P_mm;
    result_t.h = L(:);
    result_t.s = s(:);
    result_t.L = L(:);
    result_t.EVT = EVT(:);
    if ismember('date', meteo.Properties.VariableNames)
        result_t.date = meteo.date(1:Nt);
    end

    etMaxPrint = ETmax;
    if ~isscalar(etMaxPrint)
        etMaxPrint = max(etMaxPrint);
    end
    fprintf('Observed rainfall: P = %.1f mm | leaching = %.1f mm | ETmax = %.4f mm/h\n', ...
        sum(P_mm), sum(L), etMaxPrint);
    % if ismember('date', meteo.Properties.VariableNames)
    %     fprintf('Periodo meteo: %s -> %s (%d passi orari)\n', ...
    %         string(meteo.date(1)), string(meteo.date(end)), Nt);
    % end
end

function meteo = readMeteoFile(fileInput)
    if ~isfile(fileInput)
        error('generateObservedRainfall: file non trovato: %s', fileInput);
    end

    [~, ~, ext] = fileparts(fileInput);
    if any(strcmpi(ext, {'.xlsx', '.xls', '.csv'}))
        meteo = readtable(fileInput);
        if ~ismember('PRCP', meteo.Properties.VariableNames)
            error('generateObservedRainfall: tabella meteo senza colonna PRCP.');
        end
        return;
    end

    raw = importdata(fileInput);
    if isstruct(raw) && isfield(raw, 'data')
        prcp = raw.data(:);
    else
        prcp = raw(:);
    end
    meteo = table(prcp, 'VariableNames', {'PRCP'});
end

function [ETmax, ETsw] = estimateSoilET(meteo, rain, dt_h)
    etPar = rain.et;
    Kc = getRainParam(etPar, 'Kc', 0.71);
    ETsw_frac = getRainParam(etPar, 'ETsw_frac', 0.05);
    lat = getRainParam(etPar, 'latitude', 45.85);
    ETmax_default = getRainParam(etPar, 'ETmax_default', 0.5);

    cols = meteo.Properties.VariableNames;
    if ~all(ismember({'Tmax', 'Tmin'}, cols))
        ETmax = ETmax_default * Kc * 10 / 24;
        ETsw = ETsw_frac * ETmax;
        ETmax = repmat(ETmax, height(meteo), 1);
        ETsw = repmat(ETsw, height(meteo), 1);
        warning('generateObservedRainfall: Tmax/Tmin assenti; ET costante %.4f mm/h.', ETmax(1));
        return;
    end

    stepsPerDay = round(24 / dt_h);
    nDays = floor(height(meteo) / stepsPerDay);
    if nDays < 1
        ETmax = 0.5 * 10 / 24; % [mm/h] fallback (~0.5 cm/d)
        ETsw = ETsw_frac * ETmax;
        warning('generateObservedRainfall: meteo troppo corto per ET; uso ETmax = %.4f mm/h.', ETmax);
        return;
    end

    nUse = nDays * stepsPerDay;
    block = reshape(meteo.Tmax(1:nUse), stepsPerDay, nDays);
    Tmax_d = max(block, [], 1, 'omitnan')';
    block = reshape(meteo.Tmin(1:nUse), stepsPerDay, nDays);
    Tmin_d = min(block, [], 1, 'omitnan')';

    if ismember('date', meteo.Properties.VariableNames)
        doy = day(datetime(meteo.date(1:stepsPerDay:nUse), 'InputFormat', 'preserve'), 'dayofyear');
    else
        doy = (1:nDays)';
    end

    ET0_mm_d = hargreavesET0(Tmax_d, Tmin_d, lat, doy);
    ET0_mm_h = ET0_mm_d * Kc / 24;

    nMeteo = height(meteo);
    ETmax = expandDailyToHourly(ET0_mm_h, stepsPerDay, nMeteo);
    ETsw = ETsw_frac * ETmax;
end

function ET0 = hargreavesET0(Tmax, Tmin, lat_deg, doy)
    Tmean = (Tmax + Tmin) / 2;
    dT = max(Tmax - Tmin, 0);
    lat = lat_deg * pi / 180;
    dr = 1 + 0.033 * cos(2 * pi * doy / 365);
    delta = 0.409 * sin(2 * pi * doy / 365 - 1.39);
    omega_s = acos(max(min(-tan(lat) .* tan(delta), 1), -1));
    Gsc = 0.0820;
    Ra = (24 * 60 / pi) * Gsc * dr .* (omega_s .* sin(lat) .* sin(delta) + ...
        cos(lat) .* cos(delta) .* sin(omega_s));
    ET0 = 0.0023 * (Tmean + 17.8) .* sqrt(dT) .* Ra;
end

function v = expandDailyToHourly(daily, stepsPerDay, nTarget)
    v = repelem(daily(:), stepsPerDay);
    if numel(v) < nTarget
        v = [v; repmat(v(end), nTarget - numel(v), 1)];
    else
        v = v(1:nTarget);
    end
end

function soil = prepareSoilForMmH(soil)
    if isfield(soil, 'units') && strcmpi(soil.units, 'mm_h')
        return;
    end
    soil.Ks = soil.Ks * 10 / 24;
    soil.Zr = soil.Zr * 10;
end

function v = getRainParam(s, name, default)
    if isstruct(s) && isfield(s, name) && ~isempty(s.(name))
        v = s.(name);
    else
        v = default;
    end
end

function ensurePlotColors()
    global colors
    if isempty(colors) || ~isstruct(colors)
        colors = struct( ...
            'green', [0.00 0.45 0.15], ...
            'light_green', [0.38 0.85 0.30], ...
            'orange', [0.85 0.32 0.10], ...
            'light_orange', [0.93 0.46 0.20]);
    end
end
