%(proiezione nel caso 1-dim della procedura di stima della corrosione, ortogonale a S)

%Approssimazione del seguente problema:
%Dato (t,x) in [t0, tf] x [0,L], con L variabile, sia h il passo temporale
%di discretizzazione di [t0, tf] e Dx il passo di [0,L].

%equazione del calore: T=T(t,x) temperatura in [t0, tf] x [0,L]
%dT/dt - c (d^2)T/d (x^2)=0
%dT/dx(t,0)=f0(t), forzante(flash) nota

%dT/dx (t,L)=0, condizioni di Neumnann omogenee
%T(t0,x)=g0(x), condizioni iniziali note


%problema discretizzato nello spazio: Ti(t)=T(t,xi), xi=i*Dx;

% T'(t)=c A T(t)+b(t), T(t)=[T0(t),T1(t),...,Tn(t)]', b(t)=[f0(t), 0,...,0]' ; (ODE)
% T(t0)=g0 <=> [T0(t0),T1(t0),...,Tn-1(t0)]' = [g0(x0),g0(x1),...,g0(xn-1),g0(xn)]';


c=1;
tfin=10;
tin=0;
N=1000; %numero di istanti temporali
dt=(tfin-tin)/N;
save dati c tfin tin N dt;

A=sparse(n+1,n+1);     

for i=2:n %nodi interni: discretizzazione del laplaciano
    A(i,i) = -2/Dx^2;
    if i > 1      A(i,i-1) = 1/Dx^2;   end
    if i < n+1      A(i,i+1) = 1/Dx^2;   end
end

% devo imporre le CB: quindi aggiungo le equazioni corrispondenti

A(1,1)=-1/Dx; A(1,2)=1/Dx; %<=> (T1(t)-T0(t))/Dx=f0(t)
A(n+1,n+1)=-1/Dx; A(n+1,n)=1/Dx;  %CB neumann su xn <=> (Tn(t)-Tn-1(t))/-Dx=0
 
%otteniamo quindi l'equazione discretizzata nello spazio:
% T'(t)=c A T(t)+b(t), T(t)=[T0(t),T1(t),...,Tn(t)]', b(t)=[f0(t), 0,...,0]' ; (ODE)
% T(t0)=g0 <=> [T0(t0),T1(t0),...,Tn-1(t0)]' = [g0(x0),g0(x1),...,g0(xn-1),g0(xn)]';
%Osservazione: b(t) può essere pensato come l'ingresso del modello (u=b(t))

%quindi otteniamo il seguente sistema state-space:
% T'(t)=c A T(t)+b(t);
% T_vero_x0(t)=T0(t)=C T(t), C=[1,0,...,0]; <=> serie temporale delle temperature raccolte in x0

%Quindi l'ordine nx del modello state space è n+1.



%per risoluzione dell'(ODE) uso un theta metodo

T_old=g0(n,Dx); % CI (per k=0)

euler=1;    %1=Backard-Euler (EI)
            %0=Euler (EE)
save euler euler;
hist_temp=T_old;
switch euler
    case 0,
        for k=0:N-1
            b= zeros(n+1,1);
            b(1)= f0((k-1)*h+tin);
            T = dt*(c*A*T_old+b)+T_old; %<=> T(k+1)=(I+hcA)T(k)+h*b
            T_old=T;
            hist_temp=[hist_temp,T];
        end
    case 1,
        for k=0:N-1
            b= zeros(n+1,1);
            b(1)= f0(k*dt+tin);
            %we solve the system (1-dt*c*A)*T^{k+1}=dt*b+T^k
            [L,U]=lu(speye(size(A,1))-dt*c*A);
            temp = L \ (dt*b+T_old);
            T = U \ temp; %<=> T(k+1)=(I-dt*c*A)^{-1}T(k)+(I-dt*c*A)^{-1}*dt*b
            T_old=T;
            hist_temp=[hist_temp,T];
        end
end


T_vero_x0=hist_temp(1,:); % <=> T_vero_x0(k)=C T(k);


%disp('temperatura in tf nei nodi di [0,L]: [T0(tf),T1(tf),...,Tn-1(tf),Tn(tf)] = ' );
%disp(num2str(T_fin'));
