% jarkovskeho_jev_1d_regolit.m
% Jarkovskeho jev 1-D: animace teploty v hloubce; regolitovy material.
% Miroslav Broz (miroslav.broz@email.cz), Jan 28th 2006

delete(gcf);

global T_eq T_1_0_ f phi_th_ chi;

jarkovskeho_jev_params;

K = K_regolit;
rho = rho_regolit;
chi = K/(rho*C);

P = 6*3600.;		% s; denni perioda (6 h)
T_eq = 150;		% K; rovnovazna teplota
f = 1./P;
delta = sqrt(chi/(pi*f));
phi_th_ = phi_th(f,T_eq,K);
T_1_0_ = T_1_0(f,T_eq,K);

t = 0;
%t = 0.1*P;
x = 0:1e-4:3.5*delta;

T = T_x_t(x,t.*P);
h1 = image(0*t,x,transpose(T));
crange = [100 200];
caxis(crange);
set(h1, 'CDataMapping', 'scaled');
h0 = gca;
set(h0, 'XTick', []);
set(h0, 'Position', [0.175 0.08 0.575 0.82]);
colormap(hot);

hold on;
text(0, 0.5*(x(1)+x(length(x))), 'regolith', 'HorizontalAlignment', 'center');
plot([-0.5 0.5], [0 0], 'k-', 'LineWidth', 3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% stupnice

cb = colorbar;
axes(cb);
set(cb, 'Position', [0.815 0.125 0.04 0.375]);
h6 = xlabel('temperature {\itT} / K', 'FontSize', 9, 'Position', [1 92.5 1]);

axes(cb);
hold on;
tmp = T_eq;
plot([-0.5 1.5], [tmp tmp],'k:');
h2 = text(-0.5, tmp, '{\itT}_{eq}  ', 'HorizontalAlignment', 'right');

axes(h0);
hold on;
tmp = 3*delta;
plot([-0.5 0.5], [tmp tmp],'k:');
text(0.5, tmp, '  3\delta');

h5 = title([ 'time {\itt} = ' sprintf('%.1f', t/3600)  ' h  ({\itP} = ' sprintf('%.1f', P/3600) ' h)' ]);
ylabel('depth {\itx} / m');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% zarivy tok

F = T_eq + T_1_0_ * cos(2*pi*f* (t*P)); % but plot the same amplitude as T
h3 = axes;
h4 = image(F);
caxis(crange);
set(h4, 'CDataMapping','scaled');
set(h3, 'Position', [0.76 0.8 0.15 0.1]);
set(h3, 'Xtick', [], 'YTick', []);
title('flux \propto {\itF}({\itt})');
%text(1, 1, '\propto {\itF}({\itt})', 'HorizontalAlignment', 'center');

axes(h0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%break;

%% tisk

N=300;
%N=20;

set(gcf, 'PaperPosition', [0.25 2.5 4.0 6.0]);

i=0;
dt=1/N;
for t=0:1/N:1-dt+1e-8
  T = T_x_t(x,t.*P);
  F = T_eq + T_1_0_ * cos(2*pi*f* (t.*P));
  set(h1, 'CData', transpose(T));
  set(h4, 'CData', F);
  set(h5, 'String', [ 'time {\itt} = ' sprintf('%.1f', t*P/3600)  ' h  ({\itP} = ' sprintf('%.1f', P/3600) ' h)' ]);

  print('-r192','-dpng',['png/jarkovskeho_jev_1d_rez' sprintf('%05d', i) '.png'])
%  drawnow; pause(0.2);
  i=i+1
end


