function F = analyticalSolution(Lmin, Lmax, dL, tauplot, uh, beta, k, theta)
%ANALYTICALSOLUTION  f(tau,L) e F(tau,L) con IUH gamma a L = Lmax.
%
%   F = analyticalSolution(Lmin, Lmax, dL, tauplot, uh, beta, k, theta)
%
%   Condizione al bordo a Lmax:
%       f(tau, Lmax) = k^theta / Gamma(theta) * tau^(theta-1) * exp(-k*tau)
%       F(tau, Lmax) = gammainc(k*tau, theta)
%
%   Lungo le caratteristiche (stesso schema dell'esponenziale in analytical_solution):
%       F(tau,L) = gammainc(k * tau_max(tau,L), theta)
%   con tau_max da tauMaxSolver. Per theta = 1: gammainc(x,1) = 1 - exp(-x).
%
%   f(tau,L) si ottiene come dF/dtau sulla griglia (gradient), cosi' F resta in [0,1]
%   anche per theta < 1 (dove la PDF gamma ha singolarita' in tau = 0 e l'integrazione
%   diretta della f caratteristica con cumtrapz puo' superare 1).
%
%   Output: struct con campi .tau, .L, .F_L_tau, .f_L_tau (come analytical_solution).

    if nargin < 8 || isempty(theta)
        theta = 1;
    end
    if theta <= 0 || k <= 0
        error('analyticalSolution: k e theta devono essere positivi.');
    end

    tauplot = tauplot(:)';
    L = Lmin:dL:Lmax;

    NL = numel(L);
    Ntau = numel(tauplot);
    m = uh / beta;

    f_analytical = zeros(NL, Ntau);
    F_analytical = zeros(NL, Ntau);

    L_idx = unique(round(linspace(NL, 1, 4)), 'stable');
    L_idx(1) = NL;
    L_saved = L(L_idx);

    for n = 1:NL
        Li = L(n);

        for j = 1:Ntau
            tauj = tauplot(j);
            if tauj == 0
                tauj = 1e-12;
            end

            tau_max = tauMaxSolver(tauj, Li, Lmax, beta, m);
            F_analytical(n, j) = gammainc(k * tau_max, theta);
        end

        f_row = gradient(F_analytical(n, :), tauplot);
        f_row(f_row < 0) = 0;
        f_analytical(n, :) = f_row;
    end

    for n = 1:numel(L_saved)
        idx = L_idx(n);
        area_val = trapz(tauplot, f_analytical(idx, :)); %#ok<NASGU>
        % fprintf('L = %.2f m ---> area = %.8f | max F = %.6f (theta=%.3g)\n', ...
        %     L_saved(n), area_val, max(F_analytical(idx, :)), theta);
    end

    F = struct();
    F.tau = tauplot;
    F.L = L;
    F.F_L_tau = F_analytical;
    F.f_L_tau = f_analytical;
end
