%Clear all variables, globals, functions and MEX links from memory.
clear all
%closes all the open figure windows
close all
%Clear the command window and homes the cursor.
clc

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Routing methods to determine the Longitudinal Dispersion Coefficient
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% Load and plot observed data. 
% Carica i valori misurati della concentrazione nel tempo nelle sezioni 1 e
% 2
% Sezione 1=C
dataC=load('case17_section1.txt');
t_dataC=dataC(:,1)*3600;
C_dataC=dataC(:,2)/1000000;

% Sezione 2=D
dataD=load('case17_section2.txt');
t_dataD=dataD(:,1)*3600;
C_dataD=dataD(:,2)/1000000;

figure(1); 
hold on
%Plot Observed data
plot(t_dataC/3600,C_dataC,'--b','linewidth',1.5,...
        'Marker','o','MarkerSize',6,...
        'MarkerEdgeColor',[0 0 1],...
        'MarkerFaceColor',[1 1 1]);
plot(t_dataD/3600,C_dataD,'--k','linewidth',1.5,...
        'Marker','d','MarkerSize',6,...
        'MarkerEdgeColor',[1 0 1],...
        'MarkerFaceColor',[1 1 1]);
title('C(x,t)','FontSize',12,'FontWeight','bold')
% set(gca,'PlotBoxAspectRatio',[2 1 1]);
axis([min(t_dataC)/3600 max(t_dataD)/3600 0  1.3*max(C_dataC)]) 
xlabel('time [hours]','FontSize',12,'FontWeight','bold');
ylabel('c(x,t) [mg/l]','FontSize',12,'FontWeight','bold');
set(gca,'FontSize',12,'FontWeight','bold');
legend('section1','section2')
box on
grid on

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Definition and initialization of parameters and variables
% Sabine River LA
                  
M = 0.1;                % mass of dye injected, kg
Kx =9.0;                % longitudinal dispersion coefficient, m^2/s
U = 0.7;                % mean current velocity, m/s;
i_f= 17.7e-003;         % channel slope
Do = 0.85;              % mean depth, m
Bo = 13.7;              % Average width, m

% Longitudinal Coordinates of sections 1 and 2
xfix=[1130 5950]; % xcoordinate of sites along the Sabine river, m



%**************************************************************************
%**************************************************************************
%Calcola la soluzione al variare di Kx e U, cercando i valori ottimi di
%Kx e U 
%--------------------------------------------------------------------------
%Definisce la function che calcola la soluzione
fun=@simpl_CD;
%--------------------------------------------------------------------------
input_par=[Kx;U];
%--------------------------------------------------------------------------
figure(2)
hold on
% [res,fval,exitflag,output] = fminsearch(fun,input_par);
% [res,fval] = fminsearch(fun,input_par);
[res] = fminsearch(fun,input_par);
%-------------------------------------------------------------------------
% Valori ottimi calcolati di Kx e U
Kx =res(1);            % longitudinal dispersion coefficient, m^2/s
U = res(2);            % mean current velocity, m/s;
%**************************************************************************
%**************************************************************************

%%
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%Plotta la soluzione  con i valori ottimi calcolati di Kx e U
figure(3); hold on

%Calcola i tempi medi di passaggio in 1 e 2
tmean=xfix/U + 2*Kx/U^2;
t_rout=[min(t_dataC):15:max(t_dataD)]';
for tt=1:length(t_rout)
    t=t_rout(tt);
    % Calcola C(x2,t_rout(tt));
    for ttt=1:length(t_dataC)
        tau=t_dataC(ttt);
        
        %Routing a Temporal Concentration Profile with the "Hayami Solution"
        %-------------------------------------------------------------------------
        if t>tau
            Crouting(ttt)=C_dataC(ttt)*(xfix(2)-xfix(1))/(t-tau)/...
                sqrt(4*pi*Kx*(t-tau))*exp(-(xfix(2)-xfix(1)-U*...
                (t-tau))^2/(4*Kx*(t-tau)));
        else
            Crouting(ttt)=0;
        end
    end  
    CroutedD(tt)=trapz(t_dataC,Crouting);
end
CroutedDplot=CroutedD';

%Aggiunge i valori misurati
plot(t_dataC/3600,C_dataC,'--b','linewidth',1.5,...
        'Marker','o','MarkerSize',6,...
        'MarkerEdgeColor',[0 0 1],...
        'MarkerFaceColor',[1 1 1]);
plot(t_dataD/3600,C_dataD,'--b','linewidth',1.5,...
        'Marker','d','MarkerSize',6,...
        'MarkerEdgeColor',[0 0 1],...
        'MarkerFaceColor',[1 1 1]);
title('C(x,t)','FontSize',12,'FontWeight','bold')
set(gca,'PlotBoxAspectRatio',[2 1 1]);
axis([min(t_dataC)/3600 max(t_dataD)/3600 0  1.3*max(C_dataD)]) 
xlabel('time [hours]','FontSize',12,'FontWeight','bold');
ylabel('c(x,t) [mg/l]','FontSize',12,'FontWeight','bold');
set(gca,'FontSize',12,'FontWeight','bold');
box on
grid on
CroutedDplot=CroutedD';
plot(t_rout/3600,CroutedDplot,'-k','linewidth',2.5)
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%Plotta la soluzione con i valori di Kx e U forniti dal Rutherford con il 
% "Frozen Cloud Method"

%Calcola i tempi medi di passaggio in 1 e 2
tmean=xfix/U + 2*Kx/U^2;
t_rout=[min(t_dataC):15:max(t_dataD)]';
for tt=1:length(t_rout)
    t=t_rout(tt);
    %Calcola C(x2,t_rout(tt));
    for ttt=1:length(t_dataC)
        tau=t_dataC(ttt);
        %-------------------------------------------------------------------------
        %Routing a Temporal Concentration Profile with the "Frozen Cloud Method"
        Crouting(ttt)=C_dataC(ttt)*U/sqrt(4*pi*Kx*(tmean(2)-tmean(1)))*...
            exp(-U^2*((tmean(2)-tmean(1)-t+tau)^2)/...
            (4*Kx*(tmean(2)-tmean(1))));
        %-------------------------------------------------------------------------
    end  
    CroutedD(tt)=trapz(t_dataC,Crouting);
end

CroutedDplot=CroutedD';
plot(t_rout/3600,CroutedDplot,'-r','linewidth',2.5)
legend('section1','section2','Hayami','Nuovola Congelata')
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++