clear;
%Exercice 2 séance 5
DTS=load('DTS.txt');
texp=DTS(:,1); %Temps écoulé depuis l'injection du traceur (s)
Cm=DTS(:,2); %Concentration massique en traceur dans le fluide échantillonné (mg/L)
t = 0:1:3240; %Vecteur temps à haute résolution utilisé pour tracer les courbes théoriques
%%
%Calcul de tm et var avec la fonction trapz
intCm1 = trapz(texp , Cm); %Intégration sur texp de Cm
E1 = Cm./intCm1;
F1 = cumtrapz(texp , E1);
tm1 = trapz(texp , texp.*E1); %Temps de séjour moyen (s)
var1 = trapz(texp , ((texp-tm1).^2).*E1);
%Fin du calcul de tm et var avec la fonction trapz

%%
%Calcul de tm et var avec la fonction simpson
intCm2 = simpson(texp , Cm); %Intégration sur texp de Cm
E2 = Cm./intCm2;
F2 = cumtrapz(texp , E2); %Obligé de recalculer F2 car E2 != E1
tm2 = simpson(texp , texp.*E2); %Temps de séjour moyen (s)
var2 = simpson(texp , ((texp-tm2).^2).*E2);
%Fin du calcul de tm et var avec la fonction simpson

%%
%Traçage des graphes
subplot(2,2,1);
    plot(texp , E1 , 'x' , texp , E2 , 'o');
    xlabel('t (s)'), ylabel('E(t) (s^{-1})');
    legend('trapz' , 'simpson');
    title('DTS différentielle');
subplot(2,2,2);
    plot(texp , F1 , 'x' , texp , F2 , 'o');
    xlabel('t (s)'), ylabel('F(t) (s^{-1})');
    legend('trapz' , 'simpson');
    title('DTS cumulée');

%On voit que la méthode de simpson donne des résultats très proches de la
%fonction trapz. Elle est adaptée à ce type de problème.
%Fin du traçage des graphes

%%
%Calcul du tm théorique
Vr = 0.389; %Volume utile théorique du réacteur (m^3)
Qr = 55/(24*3600); %Débit nominal du réacteur (m^3.s^-1)
tmth = Vr/Qr; %Temps moyen de séjour théorique (s)

% Conclusions : On se rend compte que le temps moyen de séjour théorique
% est plus haut que le temps de séjour réel. Le débit étant fixé on peut
% donc en déduire que le volume utile du réacteur est plus faible que son
% volume utile théorique.
%Fin calcul du tm théorique

%%
%Début du NTE 1 : "Détermination d'un modèle d'association de réacteurs idéaux d'après des
%données expérimentales."

fDTS = @(n , V , texp) ( ((n*Qr/V).^n).*( (texp.^(n-1))./factorial(n-1) ).*( exp(-texp*n*Qr/V) ) );
%%Fonction anonyme prenant en entrée n et V, texp est utilisé comme vecteur
%%temps et Qr comme débit.

%Début détermination n et V optimaux et traçage des courbes
V=zeros(1,15);
var=zeros(1,15); %pré-allocation de la taille des tableaux pour accélérer l'exécution de la boucle
for n=1 : 15
    ecart = @(V) (sum((E2-fDTS(n , V , texp)).^2)); %Calcul de la variance
    [V(n) , var(n)] = fminbnd(ecart, 0 , Vr); %On récupère les différentes valeurs de V optimisées pour chaque n ainsi que l'écart induit
    %On fixe les bornes de recherche à 0 en mini car c'est le plus petit
    %volume possible et à Vr en maxi car on présume que le volume réel sera
    %légèrement inférieur.
end
[var_opt , n_opt] = min(var); %var_opt ne sera pas utilisé, n_opt est le nombre d'étage pour lequel la variance est minimale
V_opt = V(n_opt); %Valeur du volume total tel que pour n==n_opt la variance est minimale
Eth = fDTS(n_opt , V_opt , t); %Calcul de la DTS théorique avec n et V optimaux
Fth = cumtrapz(t , Eth); %Calcul de la DTS cumulée théorique
%Fin optimisation
%%

%Début traçage des courbes
subplot(2,2,3); %Traçage du E théorique et réel
    plot(texp,E1,'x' , t, fDTS(n_opt , V_opt , t))
    title(sprintf('DTS différentielle : n=%i V=%f m^{3}', n_opt, V_opt))%'DTS différentielle : n=',num2str(n_opt),' V=',num2str(V_opt),'m^{3}');
    legend('expérimental (trapz)' , 'théorique');
    xlabel('Temps (s)'), ylabel('E(t) (s^{-1})');
subplot(2,2,4); %Traçage du F théorique et réel
    plot(texp,F1,'x' , t,Fth)
    title('DTS cumulée');
    legend('expérimental (trapz)','théorique');
    xlabel('Temps (s)'), ylabel('F(t) (s^{-1})');
%Fin traçage des courbes
%Fin NTE 1

%Conclusions : Comme prévu dans les questions précédentes on trouve bien un
%volume utile plus faible que le volume théorique donné dans l'énoncé.
%%
disp 0; %Return 0 pour vérifier que le programme a été exécuté jusqu'au bout