condif_params;

gridgen;

Nx = imax;
Ny = jmax;

for ii = 1:Nx+1,
  for jj = 1:Ny+1,

    x(ii,jj) = xloc(ii);
    y(ii,jj) = yloc(jj);

    u(ii,jj) = uinlet(y(ii,jj));
    
    T(ii,jj) = Tinlet(y(ii,jj));

  end
end

% Set-up matrix to solve for convection-diffusion equation

N = (Nx+1)*(Ny+1);
A = zeros(N, N);
r = zeros(N, 1);

% First, discretize the interior points

for ii = 2:Nx,
  for jj = 2:Ny,

    kk0  = ii   + (jj-1)*(Nx+1);
    kkE  = ii+1 + (jj-1)*(Nx+1);
    kkW  = ii-1 + (jj-1)*(Nx+1);
    kkWW = ii-2 + (jj-1)*(Nx+1);
    kkN  = ii   + (jj  )*(Nx+1);
    kkS  = ii   + (jj-2)*(Nx+1);

    yp  = 0.5*(y(ii,jj+1) + y(ii,jj  ));
    ym  = 0.5*(y(ii,jj  ) + y(ii,jj-1));
    y0  = y(ii,jj);
    
    dyp = y(ii,jj+1) - y(ii,jj);
    dym = y(ii,jj  ) - y(ii,jj-1);
    dy0 = 0.5*(dyp + dym);
       
    dx = 0.5*(x(ii+1,jj) - x(ii-1,jj));

    Kp = Kinlet(yp);
    Km = Kinlet(ym);
    K0 = Kinlet(y0);
    
    if (ii == 2),
      
      r(kk0) = 0.5*u(ii,jj)/dx*(T(ii+1, jj  ) - T(ii-1, jj  )) ...
	     - K0*(T(ii+1, jj  ) - 2*T(ii  , jj  ) + T(ii-1, jj  ))/dx^2 ...
	     -( Kp*(T(ii  , jj+1) - T(ii  , jj  ))/dyp ...
             -  Km*(T(ii  , jj  ) - T(ii  , jj-1))/dym )/dy0;

      A(kk0, kk0) =                    2*K0/dx^2 + (Kp/dyp+Km/dym)/dy0;
      A(kk0, kkE) =  0.5*u(ii,jj)/dx -   K0/dx^2;
      A(kk0, kkW) = -0.5*u(ii,jj)/dx -   K0/dx^2;
      A(kk0, kkN) =                              -  Kp/dyp/dy0;
      A(kk0, kkS) =                              -         Km/dym/dy0;
      
    else

      r(kk0) = u(ii,jj)/dx*(1.5*T(ii, jj  ) - 2*T(ii-1, jj  ) + 0.5*T(ii-2, jj )) ...
	     - K0*(T(ii+1, jj  ) - 2*T(ii  , jj  ) + T(ii-1, jj  ))/dx^2 ...
	     -( Kp*(T(ii  , jj+1) - T(ii  , jj  ))/dyp ...
             -  Km*(T(ii  , jj  ) - T(ii  , jj-1))/dym )/dy0;

      A(kk0, kk0)  =  1.5*u(ii,jj)/dx + 2*K0/dx^2 + (Kp/dyp+Km/dym)/dy0;
      A(kk0, kkE)  =                  -   K0/dx^2;
      A(kk0, kkW)  = -2.0*u(ii,jj)/dx -   K0/dx^2;
      A(kk0, kkWW) =  0.5*u(ii,jj)/dx;
      A(kk0, kkN)  =                              -  Kp/dyp/dy0;
      A(kk0, kkS)  =                              -         Km/dym/dy0;
      
    end

  end
end

% Set-up Dirichlet conditions at inlet

ii = 1;
for jj = 1:Ny+1,
 
  kk0 = ii   + (jj-1)*(Nx+1);

  r(kk0     ) = T(ii,jj) - Tinlet(y(ii,jj));
  A(kk0, kk0) = 1.0;

end

% Set-up Dirichlet conditions at bottom

jj = 1;
for ii = 2:Nx+1,

  kk0 = ii   + (jj-1)*(Nx+1);

  r(kk0     ) = T(ii,jj) - Tinlet(y(ii,jj));
  A(kk0, kk0) = 1.0;

end

% Set-up Dirichlet conditions at top boundary

jj = Ny+1;
for ii = 2:Nx+1,
 
  kk0 = ii   + (jj-1)*(Nx+1);

  r(kk0     ) = T(ii,jj) - Tinlet(y(ii,jj));
  A(kk0, kk0) = 1.0;

end

% Set-up parabolized boundary conditions at outlet

ii = Nx+1;
for jj = 2:Ny,
 
  kk0 = ii   + (jj-1)*(Nx+1);
  kkW = ii-1 + (jj-1)*(Nx+1);
  kkN = ii   + (jj  )*(Nx+1);
  kkS = ii   + (jj-2)*(Nx+1);

  yp  = 0.5*(y(ii,jj+1) + y(ii,jj  ));
  ym  = 0.5*(y(ii,jj  ) + y(ii,jj-1));
  y0  = y(ii,jj);
  
  dyp = y(ii,jj+1) - y(ii,jj);
  dym = y(ii,jj  ) - y(ii,jj-1);
  dy0 = 0.5*(dyp + dym);
       
  dx = x(ii,jj) - x(ii-1,jj);

  Kp = Kinlet(yp);
  Km = Kinlet(ym);
  K0 = Kinlet(y0);
    
  r(kk0) =     u(ii,jj)/dx*(T(ii  , jj  ) - T(ii-1, jj  )) ...
         -( Kp*(T(ii  , jj+1) - T(ii  , jj  ))/dyp ...
         -  Km*(T(ii  , jj  ) - T(ii  , jj-1))/dym )/dy0;


  A(kk0, kk0) =      u(ii,jj)/dx + (Kp/dyp+Km/dym)/dy0;
  A(kk0, kkW) =     -u(ii,jj)/dx;
  A(kk0, kkN) =                  -  Kp/dyp/dy0;
  A(kk0, kkS) =                  -         Km/dym/dy0;

end

% Solve for dT

dT = -A\r;

% Update T

for ii = 1:Nx+1,
  for jj = 1:Ny+1,

    kk0 = ii   + (jj-1)*(Nx+1);
    T(ii,jj) = T(ii,jj) + dT(kk0);

  end
end

% Plot result

contour(x,y,T,21);
xlabel('x');ylabel('y');
title('Temperature contours');

% Find max wall temperature

Tmax = max(max(T(:,jlower:jupper)));

fprintf('Maximum wall temperature = %6.2f\n',Tmax);