function [ n_data, x, y, r, fxy ]  = bivariate_normal_cdf_values ( n_data )

%% BIVARIATE_NORMAL_CDF_VALUES returns some values of the bivariate normal CDF.
%
%  Discussion:
%
%    FXY is the probability that two variables A and B, which are
%    related by a bivariate normal distribution with correlation R,
%    respectively satisfy A <= X and B <= Y.
% 
%  Modified:
%
%    05 December 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    National Bureau of Standards,
%    Tables of the Bivariate Normal Distribution and Related Functions,
%    NBS, Applied Mathematics Series, Number 50, 1959.
%
%  Parameters:
%
%    Input, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.
%
%    Output, integer N_DATA, the routine increments the input value of N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, Y, the parameters of the function.
%
%    Output, real R, the correlation value.
%
%    Output, real FXY, the value of the function.
%
  n_max = 25;

  fxy_vec(1:n_max) = [ ...
    0.023, 0.155, 0.469, 0.745, 0.832, ...
    0.841, 0.138, 0.162, 0.183, 0.201, ...
    0.218, 0.234, 0.248, 0.263, 0.277, ...
    0.291, 0.305, 0.318, 0.332, 0.346, ...
    0.360, 0.374, 0.389, 0.403, 0.416 ];
  r_vec(1:n_max) = [ ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5, -0.9, -0.8, -0.7, -0.6, ...
    -0.5, -0.4, -0.3, -0.2, -0.1, ...
     0.0,  0.1,  0.2,  0.3,  0.4, ...
     0.5,  0.6,  0.7,  0.8,  0.9 ];
  x_vec(1:n_max) = [ ...
    -2.0, -1.0,  0.0,  1.0,  2.0, ...
     3.0, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2, ...
    -0.2, -0.2, -0.2, -0.2, -0.2 ];
  y_vec(1:n_max) = [ ...
     1.0,  1.0,  1.0,  1.0,  1.0, ...
     1.0,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5, ...
     0.5,  0.5,  0.5,  0.5,  0.5 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    r = 0.0;
    x = 0.0;
    y = 0.0;
    fxy = 0.0;
  else
    r = r_vec(n_data);
    x = x_vec(n_data);
    y = y_vec(n_data);
    fxy = fxy_vec(n_data);
  end