function f_obiettivo=simpl_CD(input_par)

Kx=input_par(1);
U=input_par(2);


% Definition and initialization of parameters and variables
% Bear Creek CO   
xfix=[1130 5950]; 

% Load observed data. 
% Carica i valori misurati della C(t) 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;

    %Calcola i tempi medi di passaggio in 1 e 2
    tmean = xfix/U + 2*Kx/U^2;
    %Intervallo di integrazione
    t_rout = [min(t_dataC):15:max(t_dataD)]';
    
    %Calcola l'integrale di convoluzione
    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))));
            %-------------------------------------------------------------------------
            %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
        
% Seleziona la funzione obiettivo da minimizzare    
%choice=1     %Minimizza la somma dei quadrati delle distanze tra la C calcolata e la
              %C misurata nella sezione 2
%choice=2     %Minimizza la distanza tra il picco di concentrazione
              %misurato in 2 ed il picco calcolato
choice =1;
% -------------------------------------------------------------------------
if choice==1
    %Minimizza la somma dei quadrati delle distanze tra la C calcolata e la
    %C misurata nella sezione 2
    for ttt=1:length(t_dataD)
        tau=t_dataD(ttt);
        %Calcolo il punto della griglia t_rout pił vicino a tau
        [delta_t,index]=min(abs(t_rout-tau));
        t_routCONFR(ttt)= t_rout(index);
        CroutedD_CONFR(ttt)=CroutedD(index);
    end
    %Definisce la funzione obiettivo
    
    f_obiettivo = sum((CroutedD_CONFR-C_dataD').^2)
end

% -------------------------------------------------------------------------
if choice==2
    %Individua il massimo di Concentrazione misurato
    [maxC_dataD,index_max]=max(C_dataD);
    timeC_dataD=t_dataD(index_max);
    
    %Individua il massimo di Concentrazione calcolato
    [maxCroutedD,index_max]=max(CroutedD);
    timeCroutedD=t_rout(index_max);
    
    %Definisce la funzione obiettivo: 
   
    f_obiettivo = sqrt((maxC_dataD-maxCroutedD)^2+...
        (timeC_dataD-timeCroutedD)^2)
end


hold on

plot(t_dataC/3600,C_dataC,'--b','linewidth',1.,...
        'Marker','o','MarkerSize',6,...
        'MarkerEdgeColor',[0 0 1],...
        'MarkerFaceColor',[1 1 1]);

hold on
plot(t_dataD/3600,C_dataD,'--b','linewidth',1.,...
        '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.1*max(C_dataC)]) % definisce gli assi
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,'-r','linewidth',2.)
plot(t_rout/3600,CroutedDplot,'linewidth',2.)
colormap;
legend('section1','section2','computed');

getframe;