function a = diffusivity_2d_jvb ( a0, m, omega, n, x, y )

%*****************************************************************************80
%
%% DIFFUSIVITY_2D_JVB evaluates a 2D stochastic diffusivity function.
%
%  Discussion:
%
%    The 2D stochastic diffusion equation has the form
%
%      - Del ( A(X,Y;OMEGA) Del U(X,Y;OMEGA) ) = F(X,Y).
%
%    The stochastic parameters OMEGA(1:4*M) are assumed to be independent
%    identically distributed random variables with mean value zero and 
%    variance 1.  The distribution is typically taken to be Gaussian or
%    uniform.
%
%    For physical and mathematical reasons, the diffusivity function 
%    A(X,Y;OMEGA) must be positive.  Our random KL expansion can guarantee
%    this by using an exponential format:
%
%      A(X,Y;OMEGA) = A0 + exp ( sum(j=1:4*m) OMEGA(J) * PSI(J)(X,Y) )
%
%    This function generalizes the 4-term expansion in diffusivity_2d_bnt,
%    used by Babuska, Nobile, Tempone, to a 4*M term expansion.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ivo Babuska, Fabio Nobile, Raul Tempone,
%    A stochastic collocation method for elliptic partial differential equations
%    with random input data,
%    SIAM Journal on Numerical Analysis,
%    Volume 45, Number 3, 2007, pages 1005-1034.
%
%  Parameters:
%
%    Input, real A0, the constant term in the expansion diffusion coefficient.
%
%    Input, integer M, controls the number of values OMEGA.
%    0 <= M.
%
%    Input, real OMEGA(4*M), the stochastic parameters.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), Y(N), the points where the diffusion coefficient is to 
%    be evaluated.
%
%    Output, real A(N), the value of the diffusion coefficient at (X,Y).
%

%
%  Make A the shape of X.
%  Doing this avoids problems when M=0.
%
  [ ma, na ] = size ( x );
  a = a0 * ones ( ma, na );
%
%  Compute the KL expansion of the log.
%
  arg = 0.0;

  for j = 1 : m
%
%  Reduce strength of higher order terms.
%
    jsq = 1.0 / j^2;

    arg = arg + jsq * ( ...
        omega(4*(j-1)+1) * cos ( j * pi * x ) ...
      + omega(4*(j-1)+2) * sin ( j * pi * x ) ...
      + omega(4*(j-1)+3) * cos ( j * pi * y ) ...
      + omega(4*(j-1)+4) * sin ( j * pi * y ) );

  end
%
%  Add the (positive) terms.
%
  a = a + exp ( arg );

  return
end
