% Plotting font size
fontsize=20;
% Heat capacity ratio top to bottom
Ctb=5;
% Number of heat packets
N = 80000;
% Number of steps in simulation
nstep = 15*N;
% Time scale (characteristic time)
tau=N;
%Initialize temperature-time arrays
Tt = zeros(nstep,1);
Tb = zeros(nstep,1);
% Initial temperature top
Tt(1)=1;
% Initial temperature bottom
Tb(1)=-1;
% Room temperature
Tr=-1;

for i = 2:nstep
    r = 4*rand(1,1)-2; % Random number between 2 and -2
    DT=Tt(i-1)-Tb(i-1); % temperature difference top to bottom
    if (r<DT)
        Tt(i) = Tt(i-1) - 1/N; % Move heat quanta from top to bottom
        Tb(i) = Tb(i-1) + Ctb/N;
    else
        Tt(i) = Tt(i-1) + 1/N; % Move heat quanta from bottom to top
        Tb(i) = Tb(i-1) - Ctb/N;
    end
end
figure(1)
subplot(2,1,2)
plot((1:nstep)/tau,Tt,'r',(1:nstep)/tau,Tb,'k','LineWidth',2)
xlabel('steps/N')
ylabel('2 (T-<T_0>)/\Delta T_0')
legend('Temperature top','Temperature bottom','Location','South')
ax2=gca;
ax2.FontSize=fontsize;
ax2.FontName='TimesRoman';
%%
D=load('metalblocks_lecture.asc','-ascii');
% Guessing the characteristic time
tau=100;
% Dimensionless time
t=D(:,1)/tau;
T1=D(:,2);
T2=D(:,3);
T10=mean(T1(70:90));
T20=mean(T2(70:90));
DT0=T20-T10;
Tmean0=(T20+T10)/2;
figure(1)
subplot(2,1,1)
% Plot rescaled temperature versus rescaled time
plot(t,2*(T1-Tmean0)/DT0,'k',t,2*(T2-Tmean0)/DT0,'r','LineWidth',2)
axis([0 15 -1.0 1.01])
xlabel('t/\tau')
ylabel('2 (T-<T_0>)/\Delta T_0')
legend('Temperature top','Temperature bottom','Location','South')
ax2=gca;
ax2.FontSize=fontsize;
ax2.FontName='TimesRoman';

