clear;
close all;
clc;

%--------------------------------------------------------------------------
% This code computes the global error MRE associated to both the empirical and
% the modeled mean hillslope lengths and print them on the command window.
% The intercatchments average of the MREs is also printed to cmd.
% These operations are performed over both the default interval L/sqrt(A)
% \in [0,3] and the intervals bounded by the min val inn between 3 and the
% val associated to the maximum drainage density  observed in the catchment,
% that is L/sqrt(A) \in [0, min(3,maxObserved)].
% Make sure to have run the code "p3_computeMeansAndPlot.m"
% (available at the same deposit of this code) before running this code.
%--------------------------------------------------------------------------



saveplots = false;

catchments_list_local =   ["borsoia"; "coweeta12"; "coweeta33"; "coweeta40"; "fernow14"; "fernow16"; "fernow37"; "focobon"; "montecalvello"; "poverty25"; "poverty35"; "turboloEst"; "turboloOvest"; "southForkPotts"; "rietholzback"];
catchments_maxDd_obse = [0.56; 6.07; 3.40; 5.10; 2.96; 2.42; 2.32; 2.24; 3.69; 5.11; 5.32; 7.78; 3.50; 2.00; 3.51];

maxnumToPlot = numel(catchments_list_local);

tableErr1 = nan(numel(catchments_list_local),3);
tableErr2 = nan(numel(catchments_list_local),3);

for hh_local=1:maxnumToPlot
    name_catchment = catchments_list_local(hh_local);
    disp("Processing catchment: "+upper(name_catchment));
    maxDd = catchments_maxDd_obse(hh_local);
    
    %Import and process DTM data
    orig_path = cd("Results_"+name_catchment);
    cd("Output_files");
    load("FinalResults.mat","network_Lengths","dx","totpixel","distributions_support");
    load("smoothEmpMean.mat","smoothEmpMeanKm","init_cond");
    cd(orig_path);
    xbound = maxDd*sqrt(totpixel)*dx/1000;

    % Computing analitical approximations of the mean
    % Linearized with empirical init. cond.
    konst = init_cond^2;
    handLinMean = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    handLinMean = handLinMean*dx/1000;
    % Linearized with approx. (sqrt area) init. cond.
    konst = totpixel;
    handLinMeanAppr = - network_Lengths + sqrt((network_Lengths.^2)+konst*ones(size(network_Lengths)));
    handLinMeanAppr = handLinMeanAppr*dx/1000;
    % Horton
    horton = 0.5*totpixel*dx/1000./network_Lengths;

    % Computing the error related to each analitical model
    % Specify the interval to be considered (max val for x=L/sqrt(A))
    maxval = min(xbound,3);
    % Computing errors
    errHorton = computeErr(smoothEmpMeanKm,horton,network_Lengths/sqrt(totpixel),maxval);
    errModel = computeErr(smoothEmpMeanKm,handLinMean,network_Lengths/sqrt(totpixel),maxval);
    errModelAppr = computeErr(smoothEmpMeanKm,handLinMeanAppr,network_Lengths/sqrt(totpixel),maxval);
    tableErr1(hh_local,1) = 100*errHorton;
    tableErr1(hh_local,2) = 100*errModel;
    tableErr1(hh_local,3) = 100*errModelAppr;
    % Specify the interval to be considered (max val for x=L/sqrt(A))
    maxval = 3;
    % Computing errors
    errHorton = computeErr(smoothEmpMeanKm,horton,network_Lengths/sqrt(totpixel),maxval);
    errModel = computeErr(smoothEmpMeanKm,handLinMean,network_Lengths/sqrt(totpixel),maxval);
    errModelAppr = computeErr(smoothEmpMeanKm,handLinMeanAppr,network_Lengths/sqrt(totpixel),maxval);
    tableErr2(hh_local,1) = 100*errHorton;
    tableErr2(hh_local,2) = 100*errModel;
    tableErr2(hh_local,3) = 100*errModelAppr;

end

tableErr1 = tableErr1(1:15,:);
tableErr2 = tableErr2(1:15,:);

disp("ERRORI SU L/sqrt{A} IN [0, min(3,maxObserved)]");
disp(round(tableErr1,0));
disp("ERRORI SU L/sqrt{A} IN [0, 3]");
disp(round(tableErr2,0));

catchAveraged1 = sum(tableErr1,1)/15;
catchAveraged2 = sum(tableErr2,1)/15;

disp("MEDIE INTERCATCHMENTS SU L/sqrt{A} IN [0, min(3,maxObserved)]");
disp(catchAveraged1);
disp("MEDIE INTERCATCHMENTS SU L/sqrt{A} IN [0, 3]");
disp(catchAveraged2);

disp("p4_computeMeansGlobalError.m COMPLETED his TASK.");






%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% AUXILIARY FUNCTIONS:


function err = computeErr(emp,model,xvals,maxval)

    err = 0;
    cou = 0;
    for jj=1:numel(emp)
        if (xvals(jj)<maxval)
            cou = cou + 1;
            err = err + abs((emp(jj)-model(jj))/emp(jj));
        else
            break;
        end
    end
    err = err/cou;

end