function rho = diffusivity_2d_pwc ( h, w, a, b, c, d, omega, n, x, y )

%*****************************************************************************80
%
%% DIFFUSIVITY_2D_PWC: piecewise constant diffusivity function.
%
%  Discussion:
%
%    The 2D stochastic diffusion equation has the form
%
%      - Del ( RHO(X,Y;OMEGA) Del U(X,Y;OMEGA) ) = F(X,Y).
%
%    Here, the diffusivity RHO is assumed to be a piecewise constant function,
%    defined on rectangular grid over [A,B]x[C,D] that is H cells high
%    and W cells wide.  The parameters OMEGA(H*W) are assumed to be given 
%    positive values that are bounded away from 0.
%
%    The underlying grid is assumed to be equally spaced, so each cell
%    is (D-C)/M units high and (B-A)/N wide.
%
%       ^  ^  +-------+-------+-------+-------+
%       D  H  | (H,1) | (H,2) | ----- | (H,W) |
%       |  |  +-------+-------+-------+-------+
%             | ----- | ----- | (I,J) | ----- |
%       Y  I  +-------+-------+-------+-------+
%             | (2,1) | (2,2) | ----- | (2,W) |
%       |  |  +-------+-------+-------+-------+
%       |  |  | (1,1) | (1,2) | ----- | (1,W) |
%       C  1  +-------+-------+-------+-------+
%       |  |
%       +--+--1-------------- J --------------W->
%       +--+--A-------------- X --------------B->
%
%    Indexing is tricky, since our (X,Y) coordinate system indexing and
%    conventions for indexing an array and for ordering indices do not
%    match up.  An arbitrary choice must be made, and here we associate
%    "width" and "X" and "J", as well as "height" and "Y" and "I".
%    We start in the lower left corner, and proceed right first, and
%    then move up one row.
%
%    The following coordinate systems are available:
%
%    * cell index (I,J):  (1,1) <= (I,J) <= (H,W). 
%      I refers to "height" and J to "width".
%
%    * cell counter K:    1     <= K     <= H*W
%      Count from lower left to lower right, then
%      go up one row.
%
%    * physical cell center coordinates (XC,YC)
%      the (X,Y) coordinates of the centers of a cell
%
%    * physical coordinates (X,Y) in [A,B]x[C,D].
%      the (X,Y) coordinates of any point
%      
%    * normalized coordinates (X01,Y01) in [0,1]x[0,1].
%   
%    Transformations from one coordinate system to another:
%
%    * (I,J) --> (XC,YC)
%                X01 = (2*J-1)/2/W
%                Y01 = (2*I-1)/2/H
%                XC = A + (B-A) * X01
%                YC = C + (D-C) * Y01
%
%    * (I,J) --> K
%                K = (I-1)*W+J
%
%    * K     --> (I,J)
%                I = floor ( K - 1 ) / W ) + 1
%                J = mod ( K - 1, W ) + 1
%
%    * K     --> (XC,YC)
%                I = floor ( K - 1 ) / W ) + 1
%                J = mod ( K - 1, W ) + 1
%                X01 = (2*J-1)/2/W
%                Y01 = (2*I-1)/2/H
%                XC = A + (B-A) * X01
%                YC = C + (D-C) * Y01
%
%    * (X,Y) --> (I,J)
%                X01 = ( X - A ) / ( B - A )
%                Y01 = ( Y - C ) / ( D - C )
%                I = round ( ( 2 * h * Y01 + 1 ) / 2 );
%                J = round ( ( 2 * w * X01 + 1 ) / 2 );
%
%    * (X,Y) --> K
%                X01 = ( X - A ) / ( B - A )
%                Y01 = ( Y - C ) / ( D - C )
%                I = round ( ( 2 * h * Y01 + 1 ) / 2 );
%                J = round ( ( 2 * w * X01 + 1 ) / 2 );
%                K = (I-1)*W+J
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer H, W, specifies the number of elements.
%
%    real A, B, the lower and upper limits of X.
%
%    real C, D, the lower and upper limits of Y.
%
%    real OMEGA(H*W), the parameters.
%
%    integer N, the number of evaluation points.
%
%    real X(N), Y(N), the points where the diffusion coefficient is to 
%    be evaluated.
%
%  Output:
%
%    real RHO(N), the value of the diffusion coefficient at each point.
%
  x01 = ( x - a ) / ( b - a );
  y01 = ( y - c ) / ( d - c );

  i = round ( ( 2 * h * y01 + 1 ) / 2 );
  i = max ( i, 1 );
  i = min ( i, h );

  j = round ( ( 2 * w * x01 + 1 ) / 2 );
  j = max ( j, 1 );
  j = min ( j, w );

  k = ( i - 1 ) * w + j;

  rho = omega(k);

  return
end
