getf('diffusionThermique1D.sci');
N=2000;
T=ones(N,1); 
X=linspace(0,1,N).';
for j=1:int(N*0.99),
    T(j)=2;
end;
coef=[[1,1]];
xa=0; xb=0.99;
t=0; Q=zeros(N,1); Se=zeros(N,1); Sp=zeros(N,1);
[T1,t,Q1,Se1,Sp1]=diffusion(N,"neumann",0,"neumann",0,coef,T,t,Q,Se,Sp,1e-9,1e-7);
[T2,t,Q2,Se2,Sp2]=diffusion(N,"neumann",0,"neumann",0,coef,T1,t,Q1,Se1,Sp1,1e-8,1e-6);
[T3,t,Q3,Se3,Sp3]=diffusion(N,"neumann",0,"neumann",0,coef,T2,t,Q2,Se2,Sp2,1e-7,1e-5);
[T4,t,Q4,Se4,Sp4]=diffusion(N,"neumann",0,"neumann",0,coef,T3,t,Q3,Se3,Sp3,1e-6,1e-4);
[T5,t,Q5,Se5,Sp5]=diffusion(N,"neumann",0,"neumann",0,coef,T4,t,Q4,Se4,Sp4,1e-5,1e-3);
[T6,t,Q6,Se6,Sp6]=diffusion(N,"neumann",0,"neumann",0,coef,T5,t,Q5,Se5,Sp5,1e-4,1e-2);
[T7,t,Q7,Se7,Sp7]=diffusion(N,"neumann",0,"neumann",0,coef,T5,t,Q6,Se6,Sp6,1e-3,1e-1);
             
plotA=scf();

rct=[0.9,1,1,2];
plot2d(X,T1,style=4,rect=rct)
plot2d(X,T2,style=5,rect=rct)
plot2d(X,T3,style=6,rect=rct)
plot2d(X,T4,style=9,rect=rct)
plot2d(X,T5,style=13,rect=rct)
plot2d(X,T6,style=15,rect=rct)
plot2d(X,T7,style=17,rect=rct)
xsegs([0.99,0.99],[1,2])
plotA.children.x_label.text="x";
plotA.children.y_label.text="T";
legends(["t=1e-7","t=1e-6","t=1e-5","t=1e-4","t=1e-3","t=1e-2","t=0.1"],[4,5,6,9,13,15,17])
             
plotB=scf();

rct=[0,0,1,0.01]
plot2d(X,Q7,style=5,rect=rct)            
plot2d(X,Se7,style=9,rect=rct)
plotB.children.x_label.text="x";
legends(["Q","Se"],[5,9])
             
plotC=scf();

plot2d(X,Q7,style=5,rect=[0.8,0,1,0.01])
plot2d(X,Se7,style=9,rect=[0.8,0,1,0.01])
xsegs([0.99,0.99],[0,0.01])
plotC.children.x_label.text="x";
legends(["Q","Se"],[5,9])
             
plotC1=scf();

rct=[0.8,0,1,0.01]
plot2d(X,Q1,style=4,rect=rct)
plot2d(X,Q2,style=5,rect=rct)
plot2d(X,Q3,style=6,rect=rct)
plot2d(X,Q4,style=9,rect=rct)
plot2d(X,Q5,style=13,rect=rct)
plot2d(X,Q6,style=15,rect=rct)
plot2d(X,Q7,style=17,rect=rct)
xsegs([0.99,0.99],[0,0.01])
plotC1.children.x_label.text="x";
plotC1.children.y_label.text="Q";
legends(["t=1e-7","t=1e-6","t=1e-5","t=1e-4","t=1e-3","t=1e-2","t=0.1"],[4,5,6,9,13,15,17])
             
plotD=scf();

rct=[0.95,0,1,1e-4];
plot2d(X,Sp1,style=4,rect=rct)
plot2d(X,Sp2,style=5,rect=rct)
plot2d(X,Sp3,style=6,rect=rct)
plot2d(X,Sp4,style=9,rect=rct)
plot2d(X,Sp5,style=13,rect=rct)
plot2d(X,Sp6,style=15,rect=rct)
plot2d(X,Sp7,style=15,rect=rct)
xsegs([0.99,0.99],[0,1e-4])
plotD.children.x_label.text="x";
plotD.children.y_label.text="Sp";
legends(["t=1e-7","t=1e-6","t=1e-5","t=1e-4","t=1e-3","t=1e-2","t=0.1"],[4,5,6,9,13,15,17])
             

r=zeros(N,1);
for j=int(N*0.5):int(N*0.99),
    r(j)=Q7(j)/(T7(1)*Se7(j));
end        
             
plotE=scf();

plot2d(X,r,style=5,rect=[0.5,0,1,1])
xsegs([0.99,0.99],[0,1])
plotE.children.x_label.text="x";
plotE.children.y_label.text="r";
            
j=int(0.99*N); Tm=Q7(j)/Se7(j)
