%% *****************************************
% GCH 2535 - Hiver 2019 - Devoir No2 
% Equipe:
%
%
%
%% *****************************************

clear all
% close all

%% *****************************************
%  Prétraitement
%  *****************************************

%Constantes physiques
k_eau=0.6;            %Conductivité thermique de l'eau (W/(m.K))
rho_eau=1000;         %Densité de l'eau (kg/m3)
c_eau=4186;           %Capacité thermique massique de l'eau (J/(kg.K))                 
alpha_eau=k_eau/(rho_eau*c_eau);   %Diffusivité thermique de l'eau (m2/s)
L_eau=334000;         %Entalpie de fusion de l'eau (J/kg)
h_air=5;              %Coefficient de convection thermique (W/((m2.K))
enso = 20;            %Ensoleillement (W/m2)
albd_eau = 0.02;      %Albédo des lacs (-)

%Conditions initiales
T_glace=273.15;
delta_T=10;
T_air=273.15+delta_T; %Température de l'air (K)

%Discrétisation du domaine
H=1;                  %Profondeur du domaine (m) 
h0=0;                 %Profondeur d'eau (m) 
h=0;                  %Profondeur d'eau (m) 
N=51;                 %Nombre de noeuds du domaine 
dx=H/(N-1);           %Pas en espace (m)
%Constantes de calcul liées au domaine
inv_dx=1/dx;
inv_2dx=1/(2*dx);
inv_dx_carre=inv_dx*inv_dx;

%Déclaration des vecteurs et matrices
VectB=zeros(N,1);
MatA=zeros(N,N);
VectT=zeros(N,1);
x_positions=zeros(N,1);

%Initialisation des vecteurs 
for n=1:N
    VectT(n)=T_glace;
    x_positions(n)=(n-1)*dx;
end  

%% *****************************************
%  Solveur
%  *****************************************

%Paramètres temporels
tfin=3600*24*200; % Temps final (s)
dt=2000; % pas de temps (s)
numero_iteration=0;
frequence_sortie=10;
t=0;

%Paramètres de calcul
gamma=alpha_eau*dt*inv_dx_carre;
sigma=k_eau/(rho_eau*L_eau);

%Conditions limites 
dirichlet = false;
is_source = false;

%Boucle en temps 
while (t<tfin) && (h<H)
    
    %Nombre de noeuds dans la phase liquide eau
    N_eau=floor(h/dx)+1; 
    
    %Accroissement de T par rayonnement solaire (K/s)
    if (is_source==true)
        therme_source=(1-albd_eau)*enso/(rho_eau*c_eau*h0); 
    else
        therme_source=0;
    end
    
    %Resolution du profil de température dans l'eau -----------------------
    %Conditions de Dirichlet
    if (dirichlet==true)
        MatA(1,1)=1;
        VectB(1)=T_air;        
    else
    %Conditions de Robin
        MatA(1,1)=3*k_eau*inv_2dx+h_air;
        MatA(1,2)=-4*k_eau*inv_2dx;
        MatA(1,3)=k_eau*inv_2dx;
        VectB(1)=h_air*T_air;
    end
    
    for n=2:N_eau
        MatA(n,n-1)=-gamma;
        MatA(n,n)=1+2*gamma;
        MatA(n,n+1)=-gamma;
        VectB(n)=VectT(n)+therme_source*dt;
    end   
    
    for n=N_eau+1:N
        MatA(n,n)=1;
        VectB(n)=T_glace;
    end  
    
    conditionement = cond(MatA);
    
    VectT=(MatA\VectB);
    
    %Resolution de la position de l'interface -----------------------------
    if (N_eau==0) 
        h=h0+sigma*dt*inv_dx*(T_air-T_glace);
    elseif (N_eau==1)
        h=h0+sigma*dt*inv_dx*(VectT(N_eau)-T_glace);
    else
        h=h0+sigma*dt*inv_2dx*(-VectT(N_eau-1)+4*VectT(N_eau)-3*T_glace); 
    end
    h0=h;
    
    %Sortie des résultast -------------------------------------------------
    if (mod(numero_iteration,frequence_sortie)==0)
        i = numero_iteration/frequence_sortie+1;
        h_chronologie(i)=h;
        chronologie(i)=t/3600/24;
        X = [' temps (heures) : ',num2str(t/3600),' h (m) : ',num2str(h)];
        disp(X)
    end 
    
    t=t+dt;
    numero_iteration=numero_iteration+1;
end

%% *****************************************
%  Post-traitement des résultats
%  *****************************************

% NB: Vous pouvez mettre ces graphiques à votre main et modifier cette
% partie comme il vous plaira

%Un exemple de tracage de graphique 
hFig = figure(1);
set(hFig, 'Position', [400 50 800 600])
subplot(2,1,1);
hold on
plot(x_positions,VectT,'b-','LineWidth',2);
title(['Profil de température dans l''eau après ',num2str(t/24/3600),' jours'],'fontsize',16);
xlabel('Profondeur (m)','fontsize',16);
ylabel('Temperature (K)','fontsize',16);
set(gca,'fontsize',14,'LineWidth',2)
box on 
subplot(2,1,2);
plot(chronologie,h_chronologie,'b-','LineWidth',2);
title('Position de l''interface ','fontsize',16);
xlabel('temps (jours)','fontsize',16);
ylabel('h (m)','fontsize',16);
set(gca,'fontsize',14,'LineWidth',2)
box on 

%Pour exportez vos résultats dans excel au besoin 
% outputtitle={'Profondeur (m)','Température (K)'};
% outputdata=[x_positions,vectT];
% xlswrite('PackageOutput.xls',outputtitle,'M');
% xlswrite('PackageOutput.xls',outputdata,'M','A2');