% Excercise 2
% Convective heat transfer in falling liquid films
% Analytical Methods in Heat Transfer
% Summer Semester 2013

clear
close all

Bi = 0.1;                   % Biot number
Pe = 200;                   % Peclet number
Thetaair=1;


% Definition of dimensionless coordinates in y-direction; auxiliary arrays
neta = 401;                 % number of points in y-direction
eta = linspace(0,1,neta);   % dimentionless coordinates in y-direction
deta = eta(2)-eta(1);
d1 = D1eq3_7(eta);          % matrix: discretization of the first derivative
d2 = D2eq3_8(eta);          % matrix: discretization of the second derivative
unity = eye(neta, neta);
vel = (eta -0.5*eta.*eta)'; % velocity field; weight function in orthogonality relations 

% Numerical solution of eigenvalue problem
% Eigenvalue problem: Y'' + L2*Y*(eta - 0.5*eta^2) = 0
% Discretization  in the form: A*[Y] = L2*B*[Y], where A and B are
% matrices, and [Y] is the discretized function Y
% Discretization of equation, rows from 2 to (np-1)
Amatr = d2;
Bmatr=zeros(neta,neta);
for j=2:(neta-1)
    Bmatr(j,j) = -(eta(j) - 0.5*eta(j)*eta(j));   % minus velocity vield
end
% Discretization of boundary conditions
% Condition at eta = 0, boundary condition of the first kind, first row of
% the matrix
Amatr(1,:) = unity(1,:);
% Condition at eta=1, boundary condition of the third kind, last row
Amatr(neta,:) = d1(neta,:) + Bi*unity(neta,:); 

% Eigenvalues
l2 = eig(Amatr,Bmatr);
% Eigenfunctions
[fi,void] = eig(Amatr,Bmatr);

% Elimination of artefacts and sorting
bool1 = isinf(l2);
l2corr = l2.*(1-bool1)+abs(l2).*bool1; 
[L2,il] = sort(l2corr);                     % sorting
efn = zeros(neta,neta);                         % initializing the eigenfunctions
for i=1:neta                                % sorting the eigenfunctions
    efn(:,i) = fi(:,il(i));
end

% Representation of the eigenfunctions

figure(1)
plot(eta,efn(:,1),'b*-',eta,efn(:,2),'m*-',eta,efn(:,3),'k*-',eta,efn(:,4),'g*-')
xlabel('\eta')
ylabel('First four eigenfunctions')

nmax = 100;                                  % maximal number of eigenvalues and eigenfunctions to be considered

% Orthogonality matrix, is calculated to check the orthogonality
% orthmatr = zeros(nmax,nmax);
% for i=1:nmax
%     for j=i:nmax
%         fun = efn(:,i).*efn(:,j).*vel;
%         intfun = (sum(fun)-0.5*fun(1)-0.5*fun(np))*deta; % integral, trapezoidal rule
%         orthmatr(i,j) = intfun;
%         orthmatr(j,i) = intfun;
%     end
% end
% 
% orthmatr

% Computation of coefficients


for i=1:nmax
        fun1 = efn(:,i).*efn(:,i).*vel;
        intfun1(i) = (sum(fun1)-0.5*fun1(1)-0.5*fun1(neta))*deta;
        fun2 = efn(:,i).*vel;
        intfun2(i) = (sum(fun2)-0.5*fun2(1)-0.5*fun2(neta))*deta;
        fun3 = efn(:,i).*eta'.*vel;
        intfun3(i) = (sum(fun3)-0.5*fun3(1)-0.5*fun3(neta))*deta;   
        Ccoef(i) = (-intfun3(i)*Thetaair*Bi/(Bi+1) + intfun2(i))/intfun1(i);
end

% Temperature field at certain values of ksi  - computation and graphical
% representation

ksi = [0.1 1 5 10 20 100];
nksi = length(ksi); % number of cross-sections where the temperature fild will be represented, no more than 7 
clr = 'bmgrkcy'; % line colors

figure(2)

for iksi = 1:nksi
  Theta(:, iksi) = eta'*Bi/(1+Bi)*Thetaair; % initialization

  for i=1:nmax
      Theta(:, iksi) = Theta(:, iksi) + Ccoef(i)*exp(-L2(i)*ksi(iksi)/Pe)*efn(:,i);
  end
  plot(eta, Theta(:, iksi), [clr(iksi) '-'])
  text(0.2,1.6-0.07*iksi,['ksi = ' num2str(ksi(iksi))],'Color',clr(iksi))
  hold on
end
plot(eta, eta'*Bi/(1+Bi)*Thetaair, 'k--')
text(0.5, 1.53, 'dashed line - \Theta_{\infty}')

axis([0 1 0 1.6])

xlabel('\eta')
ylabel('\Theta')

