function [ n_data, y, m, d, w ] = weekday_values ( n_data )

%*****************************************************************************80
%
%% weekday_values() returns the day of the week for various dates.
%
%  Discussion:
%
%    The CE or Common Era calendar is used, under the
%    hybrid Julian/Gregorian Calendar, with a transition from Julian
%    to Gregorian.  The day after 04 October 1582 was 15 October 1582.  
%
%    The year before 1 AD or CE is 1 BC or BCE.  In this data set,
%    years BC/BCE are indicated by a negative year value.  
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz,
%    Calendrical Calculations: The Millennium Edition,
%    Cambridge University Press, 2001,
%    ISBN: 0 521 77752 6
%    LC: CE12.R45.
%
%  Input:
%
%    integer N_DATA.  The user sets N_DATA to 0 before the first call. 
%
%  Output:
%
%    integer N_DATA.  On each call, the routine increments 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.
%
%    integer Y, M, D, the Common Era date.
%
%    integer W, the day of the week.  Sunday = 1.
%
  n_max = 34;

  d_vec = [ ...
    30, ...
     8, ...
    26, ...
     3, ...
     7, ...
    18, ...
     7, ...
    19, ...
    14, ...
    18, ...
    16, ...
     3, ...
    26, ...
    20, ...
     4, ...
    25, ...
    31, ...
     9, ...
    24, ...
    10, ...
    30, ...
    24, ...
    19, ...
     2, ...
    27, ...
    19, ...
    25, ...
    29, ...
    19, ...
     7, ...
    17, ...
    25, ...
    10, ...
    18 ];
  m_vec = [ ...
     7, ...
    12, ...
     9, ...
    10, ...
     1, ...
     5, ...
    11, ...
     4, ...
    10, ...
     5, ...
     3, ...
     3, ...
     3, ...
     4, ...
     6, ...
     1, ...
     3, ...
     9, ...
     2, ...
     6, ...
     6, ...
     7, ...
     6, ...
     8, ...
     3, ...
     4, ...
     8, ...
     9, ...
     4, ...
    10, ...
     3, ...
     2, ...
    11, ...
     7 ];
  w_vec = [ ...
    1, ...
    4, ...
    4, ...
    1, ...
    4, ...
    2, ...
    7, ...
    1, ...
    7, ...
    1, ...
    6, ...
    7, ...
    6, ...
    1, ...
    1, ...
    4, ...
    7, ...
    7, ...
    7, ...
    4, ...
    1, ...
    6, ...
    1, ...
    2, ...
    4, ...
    1, ...
    1, ...
    2, ...
    2, ...
    5, ...
    3, ...
    1, ...
    4, ...
    1 ];
  y_vec = [ ...
    - 587, ...
    - 169, ...
       70, ...
      135, ...
      470, ...
      576, ...
      694, ...
     1013, ...
     1066, ...
     1096, ...
     1190, ...
     1240, ...
     1288, ...
     1298, ...
     1391, ...
     1436, ...
     1492, ...
     1553, ...
     1560, ...
     1648, ...
     1680, ...
     1716, ...
     1768, ...
     1819, ...
     1839, ...
     1903, ...
     1929, ...
     1941, ...
     1943, ...
     1943, ...
     1992, ...
     1996, ...
     2038, ...
     2094 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    y = 0;
    m = 0;
    d = 0;
    w = 0;
  else
    y = y_vec(n_data);
    m = m_vec(n_data);
    d = d_vec(n_data);
    w = w_vec(n_data);
  end

  return
end
