function [ n_data, n, h ] = harmonic_values ( n_data )

%*****************************************************************************80
%
%% harmonic_values() returns some values of the harmonic numbers.
%
%  Discussion:
%
%    H(N) = sum ( 1 <= I < N ) 1 / I
%
%    In Mathematica:
%
%      H(N) = HarmonicNumber[n]
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 May 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N_DATA.  The user sets N_DATA to 0 before the first call.  
%    Thereafter, it should simply be the value returned by the previous 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 N, the index of a harmonic number.
%
%    real H, the value of a harmonic number.
%
  n_max = 40;

  n_vec = [ ...
     1, ...
     2, ...
     3, ...
     4, ...
     5, ...
     6, ...
     7, ...
     8, ...
     9, ...
    10, ...
    11, ...
    12, ...
    13, ...
    14, ...
    15, ...
    16, ...
    17, ...
    18, ...
    19, ...
    20, ...
    21, ...
    22, ...
    23, ...
    24, ...
    25, ...
    26, ...
    27, ...
    28, ...
    29, ...
    30, ...
    31, ...
    32, ...
    33, ...
    34, ...
    35, ...
    36, ...
    37, ...
    38, ...
    39, ...
    40 ];

  h_vec = [ ...
    1.000000000000000, ...
    1.500000000000000, ...
    1.833333333333333, ...
    2.083333333333333, ...
    2.283333333333333, ...
    2.450000000000000, ...
    2.592857142857143, ...
    2.717857142857143, ...
    2.828968253968254, ...
    2.928968253968254, ...
    3.019877344877345, ...
    3.103210678210678, ...
    3.180133755133755, ...
    3.251562326562327, ...
    3.318228993228993, ...
    3.380728993228993, ...
    3.439552522640758, ...
    3.495108078196313, ...
    3.547739657143682, ...
    3.597739657143682, ...
    3.645358704762730, ...
    3.690813250217275, ...
    3.734291511086840, ...
    3.775958177753507, ...
    3.815958177753507, ...
    3.854419716215045, ...
    3.891456753252082, ...
    3.927171038966368, ...
    3.961653797587058, ...
    3.994987130920391, ...
    4.027245195436520, ...
    4.058495195436520, ...
    4.088798225739550, ...
    4.118209990445433, ...
    4.146781419016861, ...
    4.174559196794639, ...
    4.201586223821666, ...
    4.227902013295350, ...
    4.253543038936376, ...
    4.278543038936376 ]; 

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    h = 0.0;
  else
    n = n_vec(n_data);
    h = h_vec(n_data);
  end

  return
end

