%Détermination du flux pariétal et de la température de surface
%PC2R
clear all
clc
%Données
%Nombre de pas de temps futur
r=3;
%Nombre de points d'espace
Lx=60;
%nombre de thermocouples moins la CL
ntc=2;
%Position des ntc thermocouples
%***********************
%***********************
%Capteur piston PvT LTN alpha 3
XTC(1)=0.313e-3;
XTC(2)=3.877e-3;
XTC(3)=6.1305e-3;


%Capteur piston PvT technocampus alpha 2
% XTC(1)=0.000296;
% XTC(2)=0.003784;
% XTC(3)=0.006045; 



%*******************************
%Matériau plateaux PC2R : 55NCDV7**
%*******************************
rhoc=7850*460; lambda=29;




a=lambda/rhoc;

%Lecture des données dans le fichier

don=load('C:\PC_Bureau\Mastere\Formation 2016-2017\TP Mastere 2016-2017\Essai Résine Polyester\donpiston.txt');
% [filename,pathname]=uigetfile({'*.txt','Fichier txt data  (*.txt)'});
% adresse_data=strcat(pathname,filename);
% don=load(adresse_data);
[nbl,nbc]=size(don);

%Lecture du fichier entier
temps=don(:,1);
for j=1:ntc+1,
  T(:,j)=don(:,j+1);
end

%OU BIEN 
%Diminution de la taille du fichier
% for i=1:6000, 
%    k=(i-1)*10+1;
%    temps(i)=don(k,1);
% 	for j=1:ntc+1,
%    	T(i,j)=don(k,j+1);
%   end
%   end

%Définition de la température initiale
T0=mean(T(1,:))

%Nombre de points en temps
tf=length(temps);
%Pas de temps
dt=temps(tf)/(tf-1);
%Pas d'espace
dx=XTC(ntc+1)/(Lx-1);
%Fourier
Fo=a*dt/dx^2;

%Indice correspondant à la température mesurée
for j=1:ntc,
   ind(j)=0;
   while ((dx*ind(j))<XTC(j)), ind(j)=ind(j)+1; end
   ind(j)
end

%Matrice A
A(1,1)=1+2*Fo;
A(1,2)=-2*Fo;
for j=2:Lx-1,
   A(j,j-1)=-Fo; A(j,j)=1+2*Fo; A(j,j+1)=-Fo;
end
A(Lx,Lx)=1;
AA=inv(A);

%Détermination des coefficients K
b1=2*Fo*dx/lambda;
for i=1:Lx,
   S(i)=0;
end
som=0;
for i=1:r,
   B=S; B(1)=B(1)+b1; B(Lx)=0;
   S=(AA*B')';
   for j=1:ntc,
     K(j,i)=S(ind(j))-(S(ind(j))-S(ind(j)-1))*(ind(j)-XTC(j)/dx); 
     som=som+K(j,i)*K(j,i);
  end
end
for i=1:r,
   for j=1:ntc,
      K(j,i)=K(j,i)/som;
   end
end

%Traitement du problème inverse
%Initialisation
for i=1:Lx,
   TM(i,1)=T0;
end
tn=1;
qp=0;
q(1)=0;

while (tn+r<tf)
   tps(tn)=temps(tn);
   for i=1:Lx,
      TC(i)=TM(i,tn);
   end

   b1=2*Fo*dx*qp/lambda;
   som=0;
   for k=1:r,
      B=TC; B(1)=B(1)+b1; B(Lx)=T(tn+k,ntc+1);
      TC=(AA*B')';
      for j=1:ntc,
         TJ(j)=TC(ind(j))-(TC(ind(j))-TC(ind(j)-1))*(ind(j)-XTC(j)/dx);
         e(j)=T(tn+k,j)-TJ(j);
         som=som+K(j,k)*e(j);
      end
   end
   tn=tn+1;
   q(tn)=qp+som;qp=q(tn);
   b1=2*Fo*dx*q(tn)/lambda;
   B=TM(:,tn-1); B(1)=B(1)+b1; B(Lx)=T(tn,ntc+1);
   TM(:,tn)=(AA*B);
   for j=1:ntc,
      YM(tn,j)=TM(ind(j),tn)-(TM(ind(j),tn)-TM(ind(j)-1,tn))*(ind(j)-XTC(j)/dx);
      res(tn,j)=YM(tn,j)-T(tn,j);
   end
end

tps(tn)=temps(tn);
tp=tps';
%Temps
truc (:,1)=tp;
%Flux
truc(:,2)=q';
%Température de paroi
truc(:,3)=TM(1,:)';
%Température calculée aux points de mesure
for j=4:3+ntc,
   truc(:,j)=YM(:,j-3);
end
for j=ntc+4:3+2*ntc,
   truc(:,j)=res(:,j-ntc-3);
end

%nom_de_fichier=input('nom du fichier \n','s');
%adresse_sortie=strcat(pathname,nom_de_fichier,'flux_piston.txt');
save ('C:\PC_Bureau\Mastere\Formation 2016-2017\TP Mastere 2016-2017\Essai Résine Polyester\flux_piston.txt', 'truc','-ASCII', '-TABS', '-DOUBLE')