%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Script to solve the heat equation substrate with
%%% sinusoidal T at the surface
%%% 2010 Aug 16,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Initialize variables
dz = 1e-7; %each depth step is 0.1 um
Nz = 400; % Choose the number of depth steps (should go to at least 100 um)
Nt = 5000; % Choose the number of time steps
dt = 1e-6; %Length of each time step in seconds ( 2 mS)
% the ratio of dz/dt seems to determine if we choose the appropriate
%value or not

% let Ut = K*Uzz, then K is defined as thermal conductivity/heat capacity
K =  8.4*10^-10; % "SiO2" K is kappa/C = (1.38/1.64)e-9 m^2/s

T = 0*ones(Nz+1,Nt+1);  %Create temperature matrix with Nz+1 rows, and Nt+1 columns
                         % T is 0 everywhere.
freq = 10;
period=1/freq;
time = [0:period/Nt:period]; %observe one cycle of fluctuation
 % frequency of injected power [1 Hz to 100 Hz]
T(1,:) = 1*sin(2*pi*freq*time);  %Set surface temperature

maxiter = 500;
for iter = 1:maxiter
    Tlast = T; %Save the last guess
    T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp
    for i=2:Nt+1,
        depth_2D = (T(1:end-2,i-1)-2*T(2:end-1,i-1)+T(3:end,i-1))/dz^2;
        time_1D = K*depth_2D;
        T(2:end-1,i) = time_1D*dt + T(2:end-1,i-1);
        T(end,i) = T(end-1,i); % Enforce bottom BC
    end
    err(iter) = max(abs(T(:)-Tlast(:)));  %Find difference between last two solutions
    if err(iter)<1E-4
        break; % Stop if solutions very similar, we have convergence
    end
end
if iter==maxiter;
    warning('Convergence not reached')
end

figure(1)
plot(log(err)), title('Convergence plot')
figure(2)
imagesc([0 0.1],[0 dz*Nz],T); title('Temperature plot (imagesc)')
colorbar
% figure(3)
% depth = [0:dz:Nz*dz];
% contourf(time,-depth,T); title('Temperature plot (contourf)')
% colorbar
figure(3)
%plot(time,T(1,:),time,T(21,:),time,T(41,:),time,T(61,:),time,T(81,:))
 plot(time,T(1,:),time,T(21,:),time,T(41,:),time,T(61,:),time,T(81,:))
xlabel('Time (sec)'); ylabel('Temperature (C)');
legend('0m','2um','4um','6um', '8um')