function phi = fvm_soln_SS_bdry_nodes( 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 Control 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
%%
%% My treatment of Boundary Conditions differs from Patankar
%%  I impose all boundary conditions on edges of regular control volumes
%%  This means that I don't have to use half-volumes.
%%  On the other hand, the boundary values of phi do not appear 
%%  as an explicit part of my solution 
%%  but (they can be easily recovered later)
%%
%%  THIS CODE USES PATANKAR HALF-VOLUMES
%%
%%  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 Control-volumes
%%  rhs   = right-hand-side column vector containing sources and BC terms
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% get Finite-Volume dimensions
% ------------------------------
%  NOTE: if using variable-size volumes, I need to be more careful here!
%  I could carry matrices containing 
%   dX_P, dX_E, dX_W, dx_e, dx_w, 
%   dZ_P, dZ_U, dZ_D, dz_u, dz_d, etc
  %%    [nz,nx] = size( X_P );
    %%   dx = X_edges(1,2) - X_edges(1,1);  
    %%   dz = Z_edges(2,1) - Z_edges(1,1);
  %%     dx = X_P(1,2) - X_P(1,1);  
  %%     dz = Z_P(2,1) - Z_P(1,1);
%

%  get dX, dZ - spacings between volume edges
%  -------------------------------------------
       dX = diff(X_edges')';
       dZ = diff(Z_edges);
%

%  get dX_E, dx_W, dZ_U, dZ_D - 
%          - distances to adjacent volume centers
%  -----------------------------------------------
       dX_e = diff(X_P')';
    % add dummy row at eastern edge to maintain size [nx, nz]
       dX_e = [ dX_e dX_e(:,end) ];
%
       dX_w = diff(X_P')';
    % add dummy row at western edge to maintain size [nx, nz]
       dX_w = [ dX_w(:,1) dX_w ];
%
       dZ_u = diff(Z_P);
    % add dummy row at top edge to maintain size [nx, nz]
       dZ_u = [ dZ_u(1,:)' dZ_u' ]';
%
       dZ_d = diff(Z_P);
    %  add dummy row at bottom edge to maintain size [nx, nz]
       dZ_d = [ dZ_d' dZ_d(end,:)' ]';

%

%  get f_e, f_w, f_u, f_d 
%          - fractional distances to 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 =  2 * ( K_P .* K_E ) ./( K_P + K_E );
 %%   K_w =  2 * ( K_P .* K_W ) ./( K_P + K_W );
 %%   K_u =  2 * ( K_P .* K_U ) ./( K_P + K_U );
 %%   K_d =  2 * ( K_P .* K_D ) ./( K_P + K_D );
%
    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; 
%

%
%  set Boundary Conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Specified phi
%  -------------
%   Boundary type #1 - boundaries on faces of regular finite volumes
%                      x_B = (x_U + x_P)/2
%                      -------------------
%    To eliminate phi_U from the equations
%           use   phi_u = (phi_U + phi_P)/2, 
%           i.e.  phi_U = 2 phi_u - phi_P
%    set a_U = 0 and move phi_U contribution onto b and a_P
%%  ** NOTE THAT this formulation assumes that Z_u is 
%%     at midpoint between Z_P and Z_U ** 
   %%      a_P(1,:) = a_P(1,:) + a_U(1,:); 
   %%      b(1,:) = b(1,:) + 2*phi_u .* a_U(1,:);
   %%      a_U(1,:) = 0;           % to make P independent of volume
                                   % at bottom of previous column
%
%   Boundary type #2 - boundaries are at solution points X_P (Patankar)
%    Set boundary point coefficient to unity, b to specified value
%    zero out coefficients on neighboring points
%
%
%  Specified phi gradient
%  ----------------------
%
%


%  specified flux boundaries
%  ---------------------------
%  
%   Boundary type #1 - Ed's design
%        boundaries on faces of regular finite volumes
%          e.g. on East boundary  x_B = (x_E + x_P)/2
%   ----------------------------------------------------
%     use  q_e = - K_e * (T_E - T_P)/dD_e   (flux toward E) 
%     i.e.  T_E = T_P - dx * (q_e/K_e) to eliminate T_E
     %%    a_P(:,end) = a_P(:,end) - a_E(:,end); 
%
     %%    b(:,end) = b(:,end) - dx * (q_e./K_e(:,end)) .* a_E(:,end);
     %%    a_E(:,end) = 0;
%
%   Boundary type #2 - recommended by Patankar (p. 50) 
%        boundaries are at solution points X_B = X_P
%   ------------------------------------------------
%
%
%
%   Upper boundary - specified phi = phi_u 
%   ========================================
         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  (x is positive outward)
%   =====================================================
   %  remove a_E from a_P
         a_P(:,end) = a_P(:,end) - a_E(:,end); 
%
   %  include specified boundary flux q_e in b
%         Boundary type #1 - boundaries on faces of regular finite volumes
     %%        b(:,end) = b(:,end) - dx * (q_e./K_e(:,end)) .* a_E(:,end);
%         Boundary type #2 - recommended by Patankar (p. 50) 
             b(:,end) = b(:,end) - q_e.* dZ(:,end);
%
   %  eliminate fictitious East point (outside domain)
         a_E(:,end) = 0;
%

%   West boundary - zero inward flux (x is positive inward)
%   =======================================================
   %  remove a_W from a_P
         a_P(:,1) = a_P(:,1) - a_W(:,1); 
%
   %  include specified boundary flux q_w in b
%         Boundary type #1 - boundaries on faces of regular finite volumes
    %%         b(:,1) = b(:,1) + dx * (q_w./K_w(:,1)) .* a_W(:,1);
%         Boundary type #2 - boundaries are at solution points X_P (Patankar)
             b(:,end) = b(:,end) + q_w.* dZ(:,end);
%
   %  eliminate fictitious West point (outside domain)
         a_W(:,1) = 0;
%


%   Bottom - geothermal flux
%   ========================  
   %  remove a_D from a_P
         a_P(end,:) = a_P(end,:) - a_D(end,:);
%
   %  include specified boundary flux q_d in b
%         Boundary type #1 - boundaries on faces of regular finite volumes
    %%         b(end,:) = b(end,:) - dz * (q_d./K_d(end,:)) .* a_D(end,:);
%         Boundary type #2 - boundaries are at solution points X_P (Patankar)
             b(end,:) = b(end,:) - q_d.* dX(end,:);
%
   %  eliminate fictitious Down point (outside domain)
         a_D(end,:) = 0;
%

%
%    Now call implicit solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        phi = solver( a_E, a_W, a_U, a_D, a_P, b );
%
%
%
%