function phi = fvm_soln_SS_2D_Practice_B( X_edges, Z_edges, X_P, Z_P, K_P, ...
                             S_C, S_P, phi_u, q_d, q_e, q_w )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  find steady-state phi(x,z), given heat sources S and conductivity K_P.
%  Uses Patankar Finite Volume method.
%
%  Boundary Conditions:
%  vector phi_u is fixed value of phi on upper boundary
%  vector q_w   is fixed flux on west boundary
%  vector q_e   is fixed flux on east boundary
%  vector q_d   is fixed flux on down boundary
%
% Boundary Conditions are treated using Patankar Practice B
%
%  a_P = coefficient of phi_P at point P = (i,j)
%  a_E = coefficient of phi_E at point east of (i,j)
%  a_W, a_U, a_D for points West, Up, Down
%
%  System of equations to be solved for phi
%
%       big_A * phi = rhs
%
%  big_A = matrix with coefficients a_E etc, modified to include BC
%  phi   = unknown values in finite volumes
%  rhs   = right-hand-side column vector containing sources and BC terms
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% get Finite-Volume dimensions
% ------------------------------
%  get dX, dZ - spacings between volume edges
%    diff function takes differences between adjacent rows, then
%    need to drop 1 row or column to get numbers of volumes 
%  --------------------------------------------------------
       dX = diff(X_edges, 1, 2);
       dZ = diff(Z_edges, 1, 1);
%
%  get dX_e, dx_w, dZ_u, dZ_d - 
%          - distances to adjacent volume centers
%  -----------------------------------------------
       dX_e = diff(X_P, 1, 2);   
    % add dummy row at eastern edge to maintain correct number
    %   of volumes [nx+1, nz+1]
       dX_e = [ dX_e dX_e(:,end) ];
%
       dX_w = diff(X_P, 1, 2);
    % add dummy row at western edge to maintain correct number
    %   of volumes [nx+1, nz+1]
       dX_w = [ dX_w(:,1) dX_w ];
%
       dZ_u = diff(Z_P, 1, 1);
    % add dummy row at top edge to maintain correct number
    %   of volumes [nx+1, nz+1]
       dZ_u = [ dZ_u(1,:)' dZ_u' ]';
%
       dZ_d = diff(Z_P, 1, 1);
    %  add dummy row at bottom edge to maintain correct number
    %   of volumes [nx+1, nz+1]
       dZ_d = [ dZ_d' dZ_d(end,:)' ]';
%
%
%  get f_e, f_w, f_u, f_d 
%          - fractional distances to centers of volume interfaces
%  --------------------------------------------------------------
       f_e = ones(size(X_P)) - (X_edges(:,2:end) - X_P)./ dX_e;
       f_w = ones(size(X_P)) - (X_P - X_edges(:,1:(end-1)) )./ dX_w;
       f_u = ones(size(Z_P)) - (Z_P - Z_edges(1:(end-1),: ))./ dZ_u;
       f_d = ones(size(Z_P)) - (Z_edges(2:end,: ) - Z_P)./ dZ_d;


% set up b matrix - includes constant part of source matrix S 
%         S = S_C + S_P * phi_P    (units are W m^-3)
% ----------------------------------------------------------
         b = S_C.* dX.* dZ;
%
%
% set up matrices for conductivities in adjacent volumes
% ------------------------------------------------------
    K_E   =  [ K_P(:,2:end) K_P(:, end) ];
    K_W   =  [ K_P(:,1) K_P(:,1:end-1) ];
    K_U   =  [ K_P(1,:)' K_P(1:end-1,:)' ]';
    K_D   =  [ K_P(2:end,:)' K_P(end,:)' ]';
%
%
% set up interface conductivity matrices
%   (These expressions need to be modified when dx and dz
%    are not uniform)
% --------------------------------------
    K_e =  1./ ( (1 - f_e)./K_P + f_e./K_E );
    K_w =  1./ ( (1 - f_w)./K_P + f_w./K_W );
    K_u =  1./ ( (1 - f_u)./K_P + f_u./K_U );
    K_d =  1./ ( (1 - f_d)./K_P + f_d./K_D );
%
%
% set up coefficients at adjacent points
% ---------------------------------------
    a_E = (K_e ./ dX_e).* dZ;
    a_W = (K_w ./ dX_w).* dZ;
    a_U = (K_u ./ dZ_u).* dX;
    a_D = (K_d ./ dZ_d).* dX;

         
% set up coefficient at central point P (incl source dependent on phi)
% --------------------------------------------------------------------
    a_P = a_E + a_W + a_U + a_D - S_P.* dX.* dZ; 
%
%
%  Include Boundary Conditions
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Upper boundary - specified phi = phi_u 
% -----------------------------------------
%    Set boundary point coefficient to unity, b to specified value
%    zero out coefficients on neighboring points
%   ========================================
         a_P(1,:) = 1; 
         b(1,:) = phi_u;
         a_U(1,:) = 0; 
         a_E(1,:) = 0;
         a_W(1,:) = 0;
         a_D(1,:) = 0;
%
%
%   East boundary - outward flux q_b (x is positive outward)
%   ========================================================
%         K_w ( phi_B - phi_W )/dX_B = q_e 
%  Set a_P = a_W
%  set b   = -q_e.*dZ  (dZ needed because a_W includes dZ)
% --------------------
        a_P(:,end) = a_W(:,end); 
   %
   %  zero out other coefficients
         a_E(:,end) = 0;
         a_U(:,end) = 0;
         a_D(:,end) = 0;
   % put in boundary flux (integrated over depth)
         b(:,end) = -q_e.*dZ(:,end);
%
%
%   West boundary - zero inward flux (x is positive inward)
%   =======================================================
%          K_e ( phi_E - phi_B )/dX_B = q_w 
%   So set a_P = a_E
%      set b   = q_b.*dZ  (dZ needed because a_E includes dZ)
% --------------------
        a_P(:,1) = a_E(:,1); 
   %
   %  zero out other coefficients
         a_W(:,1) = 0;
         a_U(:,1) = 0;
         a_D(:,1) = 0;
   % put in boundary flux (integrated over depth)
         b(:,1) =  q_w.*dZ(:,1);
   %
%   Bottom - geothermal flux
%   ========================  
%          K_d ( phi_B - phi_U )/dZ_B = q_d 
%   So set a_P = a_U
%      set b   = -q_d.*dX  (dX needed because a_U includes dX)
% --------------------
        a_P(end, :) = a_U(end, :); 
   %
   %  zero out other coefficients
         a_W(end, :) = 0;
         a_E(end, :) = 0;
         a_D(end, :) = 0;
   % put in boundary flux (integrated over x extent of volume)
         b(end, :) =  -q_d.*dX(end,:);

%
%  Matrix will be singular because the corner points have all zero
%  coefficients. So, set corner values to average of 2 neighbors.
%  (We could get fancier if node separations are unequal.)
%  (Don't forget to also zero out the source vector. :-)
%  ---------------------------------------------------------------
%
%  test whether to set corners
   %%  set_corners = 0;   % don't set corner equations
     set_corners = 1;   % set corner equations
     
     if(set_corners == 1)
%
 %    upper left  
        a_P(1,1) = 1;
        a_E(1,1) = 1/2;
        a_D(1,1) = 1/2;
        b(1,1)   = 0;
 %
 %    upper right
        a_P(1,end) = 1;
        a_W(1,end) = 1/2;
        a_D(1,end) = 1/2;
        b(1,end)   = 0;
 %
 %    lower left
        a_P(end,1) = 1;
        a_E(end,1) = 1/2;
        a_U(end,1) = 1/2;
        b(end,1)   = 0;
 %
 %    lower right
        a_P(end,end) = 1;
        a_W(end,end) = 1/2;
        a_U(end,end) = 1/2;
        b(end,end)   = 0;
%
     end   %  if(set_corners == 1) ...
%
%   Now call implicit solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        phi = solver( a_E, a_W, a_U, a_D, a_P, b );
%
%
%
%