% Excercise 1
% 1D transient heat conduction in a plate
% Analytical Methods in Heat Transfer
% Summer Semester 2013

clear
close all

Bi = 0.1;           % Biot number

% Determination of eigenvalues
% Graphical representation of equation, from which the eigenvalues can be
% determined
% The points of intersection between the green and blue curves correspond
% to eigenvalues
lammax = 5*pi; % the maximal value of lambda to be represented in the plot
Lam = linspace(pi/400, lammax, 100);
figure(1)
plot(Lam, Bi./Lam, '-g', Lam, tan(Lam), '-b')
hold on
axis([0 lammax -3 3]) % 
grid
xlabel('\Lambda')

% Choose between the graphical method (using ginput function) and numerical
% method (using fzero function)

% % Graphical method
% matr=ginput;
% L = matr(:,1);    % L's are the eigenvalues
% nL = length(L);

% Numerical method (function fzero)
% The function fzero is used in the form, in which the last argument is the
% initial guess of the solution. You may want to use another form, where
% the last argument is a vector containing two values. These values bound
% an interval within which the solution is sought for.
nL = 15; % the number of eigenvalues we are looking for
L = zeros(nL,1);
L(1) = fzero(@ (z) tan(z)-Bi/z, sqrt(Bi));
for i=2:nL
    L(i) = fzero(@ (z) tan(z)-Bi/z, pi*(i-1)); % L's are the eigenvalues
end

figure(1)
plot(L, Bi./L, 'r*')


figure(2) 
nZ = 50;
Z = linspace(0,1,nZ); % dimensionless coordinate
Fo = [0 0.01 0.1 0.25 0.5 1 5 10]; % this array contains the Fourier numbers (dimensionless time)
nFo = length(Fo);
tmp = ones(nZ, nFo); % initializing the temperature
for iFo=1:nFo
  tmp(:,iFo) = zeros(nZ,1);
  for iL=1:nL
      C(iL) = 2*sin(L(iL))/L(iL)/(1+(sin(L(iL)))^2/Bi);
    term = C(iL)*exp(-L(iL)^2*Fo(iFo))*cos(L(iL)*Z);
    tmp(:,iFo) = tmp(:,iFo) + term';
  end
  plot(Z,tmp(:,iFo))
  hold on
end

xlabel('\xi')
ylabel('\theta')
% Each curve corresponds to a single time instant, or a value of Fo
