function phi = fdm_soln_SS( 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 = S_C + S_P phi.
%%  Uses Finite Difference method.
%%
%%  Boundary Conditions:
%%  vector phi_u is fixed value of phi on upper boundary
%%  vector q_w   is fixed flux on west bou
%%  vector q_e   is fixed flux on east boundary
%%  vector q_d   is fixed flux on down boundary
%%
%%  d/dx[ K d phi/dx] + d/dz[ K d phi/dz] + S = 0 
%%
%%  or
%%
%%  K ( d^2phi/dx^2 + d^2phi/dz^2 ) + ...
%%           ( dK/dx dphi/dx + dK/dz dphi/dz ) + S = 0 
%%
%%  All derivatives are estimated at the grid points, using
%%  values at adjacent grid points.
%%  There are no finite boxes to consider.
%%
%%
%%  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
%%
%%
%%  a_P phi_P = a_E phi_E + a_W phi_W + a_U phi_U + a_D phi_D + b
%%
%% or:
%%
%%  - a_P phi_P + a_E phi_E + a_W phi_W + a_U phi_U + a_D phi_D = - b
%%
%%  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 at grid points (X_P, Z_P)
%%  rhs   = right-hand-side column vector - sources "b" and BC terms
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% get grid dimensions - we assume uniform dx, dz spacing
% -------------------------------------------------------
       [nz,nx] = size( X_P );
 
       dx = 1;      % for 1
       if(nx > 1 )
          dx = ( X_edges(1,end) - X_edges(1,1) )/(nx-1);
       end   %  if(nx > 1 ) ...
       
       dz = Z_edges(2,1) - Z_edges(1,1);
%  

% set up RHS vector - includes constant part of source matrix S 
%         S = S_C + S_P * phi_P    (units are W m^-3)
% ----------------------------------------------------------
         rhs = S_C;
%
%
% 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 conductivity gradient matrices
% --------------------------------------
    dK_dx =  (K_E - K_W ) ./(2*dx);
    dK_dz =  (K_D - K_U ) ./(2*dz);
%

% set up coefficients at adjacent points
% ---------------------------------------

% 1-D 
   a_E = zeros( size(K_P) );
   a_W = zeros( size(K_P) );
   
   if(nx > 1)
      a_E = ( K_P./dx.^2 + dK_dx./(2*dx) );
      a_W = ( K_P./dx.^2 - dK_dx./(2*dx) );
   end
    a_D = ( K_P./dz.^2 + dK_dz./(2*dz) );
    a_U = ( K_P./dz.^2 - dK_dz./(2*dz) );

         
% set up coefficient at central point P (incl source dependent on phi)
% --------------------------------------------------------------------
  %%  a_P = 2*K_P.*(1./dx.^2 + 1./dz.^2) - S_P;
    a_P = a_E + a_W + a_U + a_D - S_P; 
%
%

%  set Boundary Conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    specified phi = phi_u on upper boundary
%    ----------------------------------------
%           use   phi_P = phi_u; 
%    i.e. set a_P = 1, 
%         set a_U = a_D = a_E = a_W = 0
%         set rhs = phi_u
%%  ** NOTE THAT this FD formulation doesn't actually ever use 
%%  ** source term at points on boundary where phi is specified
%
         a_P(1,:) = 1; 
         a_U(1,:) = 0;
         a_D(1,:) = 0;
         a_E(1,:) = 0;
         a_W(1,:) = 0;
         rhs(1,:) = phi_u;
%

%    specified flux boundaries
%    ---------------------------
%     East boundary: q_e = outward flux 
%           use  q_e = - K_P * (T_E - T_W)/(2 dx)   (flux toward E) 
%          i.e.  T_E = T_W - 2*dx * (q_e/K_P)       to eliminate T_E 
%
         rhs(:,end) = rhs(:,end) - 2*dx * a_E(:,end).* (q_e./K_P(:,end) );
         a_W(:,end) = a_W(:,end) + a_E(:,end);
         a_E(:,end)   = 0;
%        
  
         
%     West boundary: q_w = inward flux 
%           use  q_w = -K_P * (T_E - T_W)/(2 dx)  (flux from W)
%          i.e.  T_W = T_P + 2*dx * (q_w/K_P)      to eliminate T_W
%
         rhs(:,1) = rhs(:,1) + 2*dx * a_W(:,1).* (q_w./K_P(:,1));
         a_E(:,1) = a_E(:,1) + a_W(:,1);
         a_W(:,1)   = 0;

%     Bottom boundary: q_d = geothermal flux (upward, inward)
%           use  q_d = -K_P * (T_D - T_U)/(2 dz)    (flux down)
%          i.e.  T_D = T_U - 2*dz * (q_d / K_P)      to eliminate T_D
%            must also modify rhs
         rhs(end,:) = rhs(end,:) - 2*dz * a_D(end,:).* (q_d./K_P(end,:));
         a_U(end,:) = a_U(end,:) + a_D(end,:);
         a_D(end,:) = 0;
%    
%
%    Now call matrix solver
         phi = solver( a_E, a_W, a_U, a_D, a_P, rhs );
%
%
%
%