% Octave script m file
% solve one dimension  diffusion partial differential eq. by explicit method
% main script
% Constants for Calculation
D=1.75; % cm^2/s diffusion constant
length=10; % cm
N=100; % Divsion Number of the length
lambda=1/6;
os2lambda=1-2*lambda;
TR=30000; % Time Range
Tint=100; % Time interval
Dint=Tint*10; % Time interval for final graph
SaveFname="Tsave.txt"; % File name to save data
fprintf("Simulation of Temperature distribution in one dimensional thermal diffusion\n");
fprintf("  Thermal Diffusion Constant : %7.4f cm^2/s \n", D);
fprintf("  Length : %7.4f cm\n",length);
fprintf("  Number of divided of x range : %d \n",N);
fprintf("  Value of parameter lambda : %7.4f \n",lambda);
fprintf("  Number of time iteration : %d \n", TR);
fprintf("  File name to save data : %s\n", SaveFname);
%
N1=N+1;
TR1=TR+1;
%
%
DeltaX=length/N;
Deltat=lambda*DeltaX^2/D;
X=[0:DeltaX:length]';
%
% construct a sparse Matrix
dia=[1 (ones(1,N-1)*os2lambda) 1];
Cd=sparse(1:N+1,1:N+1,dia,N+1,N+1);
Cu=sparse(2:N,3:N+1,lambda,N+1,N+1);
Cl=sparse(2:N,1:N-1,lambda,N+1,N+1);
C=Cd+Cu+Cl;
% Initial and boundary condition
Tp=[30 (ones(1,N-1)*50) 30]';  %Initial condition
Tb1=ones(TR1+1,1)*30;
Tb2=Tb1;
%
Tsave = fopen (SaveFname, "w");
% solve PDE by iteration
for m=1:TR1
  if (mod((m-1),Tint) == 0) 
	xlabel("x / cm\n");  ylabel("T / deg\n"); axis([0 10 30 50]);
	plot(X,Tp, sprintf(";%d: t=%7.4f s;",(m-1),(m-1)*Deltat));
	xlabel("x / cm\n");  ylabel("T / deg\n"); axis([0 10 30 50]);
    drawnow;  % required for Octave 2.9.xxx
    fprintf(Tsave, "%e ",Tp); fputs(Tsave,"\n"); 
  end
  Tn=C*Tp;
  Tp=Tn;
end
fclose(Tsave);
input("Press any key to continue");
Tread = fopen (SaveFname, "r");
Tx=zeros(1,N1);
clf;
xlabel("x / cm\n");  ylabel("T / deg\n");  axis([0 10 30 50]);
hold on
for m=1:Tint:TR1
   sTx=[ "[" fgets(Tread) "]" ];
   Tx=eval(sTx);
  if (mod((m-1),Dint) == 0)
    clr=mod(floor((m-1)/Dint),5)+1;
	;
    plot(X,Tx', sprintf("%d;%d: t=%7.4f s;",clr,(m-1),(m-1)*Deltat));
  end
end
fclose(Tread);