function [ n_data, jed, mjd ] = jed_mjd_values ( n_data )

%*****************************************************************************80
%
%% jed_mjd_values() returns the MJD for Julian Ephemeris Dates.
%
%  Discussion:
%
%    The JED (Julian Ephemeris Date) is a calendrical system which counts days,
%    starting from noon on 1 January 4713 BCE.
%
%    The MJD (Modified Julian Day) counts days starting from midnight, 
%    17 November 1858.  This essentially subtracts 2400000.5 days from the JED.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold and Nachum Dershowitz,
%    Calendrical Calculations: The Millennium Edition,
%    Cambridge University Press, 2001,
%    ISBN: 0 521 77752 6
%
%  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.
%
%    real JED, the Julian Ephemeris Date.
%
%    real MJD, the Modified Julian Ephemeris Date.
%
  n_max = 33;

  jed_vec = [ ...
     1507231.5E+00, ...
     1660037.5E+00, ...
     1746893.5E+00, ...
     1770641.5E+00, ...
     1892731.5E+00, ...
     1931579.5E+00, ...
     1974851.5E+00, ...
     2091164.5E+00, ...
     2121509.5E+00, ...
     2155779.5E+00, ...
     2174029.5E+00, ...
     2191584.5E+00, ...
     2195261.5E+00, ...
     2229274.5E+00, ...
     2245580.5E+00, ...
     2266100.5E+00, ...
     2288542.5E+00, ...
     2290901.5E+00, ...
     2323140.5E+00, ...
     2334848.5E+00, ...
     2348020.5E+00, ...
     2366978.5E+00, ...
     2385648.5E+00, ...
     2392825.5E+00, ...
     2416223.5E+00, ...
     2425848.5E+00, ...
     2430266.5E+00, ...
     2430833.5E+00, ...
     2431004.5E+00, ...
     2448698.5E+00, ...
     2450138.5E+00, ...
     2465737.5E+00, ...
     2486076.5E+00 ];

  mjd_vec = [ ...
     -892769.0E+00, ...
     -739963.0E+00, ...
     -653107.0E+00, ...
     -629359.0E+00, ...
     -507269.0E+00, ...
     -468421.0E+00, ...
     -425149.0E+00, ...
     -308836.0E+00, ...
     -278491.0E+00, ...
     -244221.0E+00, ...
     -225971.0E+00, ...
     -208416.0E+00, ...
     -204739.0E+00, ...
     -170726.0E+00, ...
     -154420.0E+00, ...
     -133900.0E+00, ...
     -111458.0E+00, ...
     -109099.0E+00, ...
      -76860.0E+00, ...
      -65152.0E+00, ...
      -51980.0E+00, ...
      -33022.0E+00, ...
      -14352.0E+00, ...
       -7175.0E+00, ...
       16223.0E+00, ...
       25848.0E+00, ...
       30266.0E+00, ...
       30833.0E+00, ...
       31004.0E+00, ...
       48698.0E+00, ...
       50138.0E+00, ...
       65737.0E+00, ...
       86076.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    jed = 0.0;
    mjd = 0.0;
  else
    jed = jed_vec(n_data);
    mjd = mjd_vec(n_data);
  end

  return
end
