% use Matlab PDE solver to sove 1-D heat conduction
% some dimensions
k = 60.5;          % W/m-K  low carbon steel
alpha = 17.7e-6;   % m^2/s   low carbon steel
q0 = 10000;          % W/m^2
L_d = 0.050;       % m
dt_d = 8.5;        % s
dt = dt_d*alpha/L_d^2;
t = [ -.2407 : dt : 2.4 ]';

L = 1;
x = linspace(0,L,20);
% t = [linspace(0,0.5,20), linspace(0.65,1.5,9)];
% t = [ -0.24:0.06:1.8 ]';
% t = [ -0.24:0.06:2.40 ]';
m = 0;  % 0=cartesian; 1=radial;  2=spherical
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
colormap hot
% pcolor(x,t,sol)
surf(x,t,sol)
colorbar
xlabel('Distance x','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')
figure(97); clf; 
plot(t, sol(:,1), t, sol(:,end) )
legend('x=0', 'x=L' )
% figure(98); clf;
% qbc = zeros(length(t),1);
% for it = 1:length(t)
%     qbc(it) = q_of_t( t(it) );
% end
% plot( t, qbc );
% T_1pulse( t, 0, 1.2, x(1) );
% figure(97)
% % plot( t, sol(:,1) - 0.6*T_1pulse( t, 0, 1.2, x(1) )', 'ro-' )
% plot( t, sol(:,1), 'ro-' )
xlabel('time'); ylabel('temperature');
figure(99); clf;
plot(t, qspline(t), 'ro-' );
xlabel('time'); ylabel('heat flux')
TL = sol(:,end);
qbc = qspline(t);
save HeatEqOut t TL qbc
TL_d = TL*q0*L_d/k;
figure(89); clf;
plot(t*L_d^2/alpha, qspline(t)*q0/1000, 'ro-' );
xlabel('time, seconds'), ylabel('Heat Flux, kW/m^2');
figure(87); clf; 
plot(t*L_d^2/alpha, TL_d,'ro-' )
xlabel('time, seconds'), ylabel('Temperature rise, C');
t_d = t*L_d^2/alpha;
qbc_d = qbc * q0;
save HeatEqOut_d t_d TL_d qbc_d
% print a tableau
fprintf('%5.2f %5.2f %6.1f\n',[ t_d TL_d qbc_d ]' )






function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;      % coefficient of du/dt 
f = dudx;   % any combination of x,t,u,du/dx; e.g. (1.2*u + 0.3)*dudx
s = 0;      % source - any combination of x,t,T,du/dx
end
function u0 = heatic(x)
u0 = 0.0;   % any function of x
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
% p(x,t,u) + q(x,t)*f(x,t,u,dudx) = 0
%pl = q_of_t(t);
pl = qspline(t);
ql = 1;
pr = 0;
qr = 1;
end
function [q] = q_of_t( t )
% define the boundary condition at x=0
if(t <= 0)
    q = 0;
elseif( t <= 0.6 )
    q = t;
elseif( t<= 1.2 )
    q = (1.2 - t);
else
    q = 0;
end
end
function [qfine] = qspline(tfine)
tbase = [ -0.25 0. 0.15  .25 .5 .75  1. 1.25 1.5 1.75  2.0 2.1  2.2 2.25 2.3 2.4 2.5]';
vals =  [    0. 0. 0.1  .3 .7 1.0  0.7 0.6 0.4  0.2  0.1 0.02 0.01 0.0 0.0 0.0 0.0]';
% qfine = spline(tbase, vals, tfine);
qfine = pchip(tbase, vals, tfine);
end