function [ f, g ] = r89sifg ( x )

%*****************************************************************************80
%
%% r89sifg() evaluate terms needed in an approximation to SI(X) for 4 <= |X|.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    12 January 2021
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    FORTRAN90 version by John Burkardt.
%
%  Input:
%
%    real X, the point at which the terms are needed.
%    |X| must be at least 4.
%
%  Output:
%
%    real F, G, the values of F and G at X.
%
  f1cs = [ ...
  -0.1191081969051363610, -0.0247823144996236248, ...
   0.0011910281453357821, -0.0000927027714388562, ...
   0.0000093373141568271, -0.0000011058287820557, ...
   0.0000001464772071460, -0.0000000210694496288, ...
   0.0000000032293492367, -0.0000000005206529618, ...
   0.0000000000874878885, -0.0000000000152176187, ...
   0.0000000000027257192, -0.0000000000005007053, ...
   0.0000000000000940241, -0.0000000000000180014, ...
   0.0000000000000035063, -0.0000000000000006935, ...
   0.0000000000000001391, -0.0000000000000000282 ];
%
% series for f2   on the interval  0.00000 to  2.00000D-02
%                                    with weighted error   4.32D-17
%                                     log weighted error  16.36
%                           significant figures required  14.75
%                                decimal places required  17.10
%
  f2cs = [ ...
    -0.0348409253897013234, ...
    -0.0166842205677959686, ...
     0.0006752901241237738, ...
    -0.0000535066622544701, ...
     0.0000062693421779007, ...
    -0.0000009526638801991, ...
     0.0000001745629224251, ...
    -0.0000000368795403065, ...
     0.0000000087202677705, ...
    -0.0000000022601970392, ...
     0.0000000006324624977, ...
    -0.0000000001888911889, ...
     0.0000000000596774674, ...
    -0.0000000000198044313, ...
     0.0000000000068641396, ...
    -0.0000000000024731020, ...
     0.0000000000009226360, ...
    -0.0000000000003552364, ...
     0.0000000000001407606, ...
    -0.0000000000000572623, ...
     0.0000000000000238654, ...
    -0.0000000000000101714, ...
     0.0000000000000044259, ...
    -0.0000000000000019634, ...
     0.0000000000000008868, ...
    -0.0000000000000004074, ...
     0.0000000000000001901, ...
    -0.0000000000000000900, ...
     0.0000000000000000432 ];
%
% series for g1   on the interval  2.00000D-02 to  6.25000D-02
%                                    with weighted error   5.48D-17
%                                     log weighted error  16.26
%                           significant figures required  15.47
%                                decimal places required  16.92
%
  g1cs = [...
    -0.3040578798253495954, ...
    -0.0566890984597120588, ...
     0.0039046158173275644, ...
    -0.0003746075959202261, ...
     0.0000435431556559844, ...
    -0.0000057417294453025, ...
     0.0000008282552104503, ...
    -0.0000001278245892595, ...
     0.0000000207978352949, ...
    -0.0000000035313205922, ...
     0.0000000006210824236, ...
    -0.0000000001125215474, ...
     0.0000000000209088918, ...
    -0.0000000000039715832, ...
     0.0000000000007690431, ...
    -0.0000000000001514697, ...
     0.0000000000000302892, ...
    -0.0000000000000061400, ...
     0.0000000000000012601, ...
    -0.0000000000000002615, ...
     0.0000000000000000548 ];
%
% series for g2   on the interval  0.00000 to  2.00000D-02
%                                    with weighted error   5.01D-17
%                                     log weighted error  16.30
%                           significant figures required  15.12
%                                decimal places required  17.07
%
  g2cs = [
    -0.0967329367532432218, ...
    -0.0452077907957459871, ...
     0.0028190005352706523, ...
    -0.0002899167740759160, ...
     0.0000407444664601121, ...
    -0.0000071056382192354, ...
     0.0000014534723163019, ...
    -0.0000003364116512503, ...
     0.0000000859774367886, ...
    -0.0000000238437656302, ...
     0.0000000070831906340, ...
    -0.0000000022318068154, ...
     0.0000000007401087359, ...
    -0.0000000002567171162, ...
     0.0000000000926707021, ...
    -0.0000000000346693311, ...
     0.0000000000133950573, ...
    -0.0000000000053290754, ...
     0.0000000000021775312, ...
    -0.0000000000009118621, ...
     0.0000000000003905864, ...
    -0.0000000000001708459, ...
     0.0000000000000762015, ...
    -0.0000000000000346151, ...
     0.0000000000000159996, ...
    -0.0000000000000075213, ...
     0.0000000000000035970, ...
    -0.0000000000000017530, ...
     0.0000000000000008738, ...
    -0.0000000000000004487, ...
     0.0000000000000002397, ...
    -0.0000000000000001347, ...
     0.0000000000000000801, ...
    -0.0000000000000000501 ];

  eta = 0.1 * eps;
  nf1 = inits ( f1cs, 20, eta );
  nf2 = inits ( f2cs, 29, eta );
  ng1 = inits ( g1cs, 21, eta );
  ng2 = inits ( g2cs, 34, eta );

  xbig = sqrt ( 1.0 / eps );
  xmaxf = exp ( min ( - log ( realmin ( ) ), log ( realmax ( ) ) ) - 0.01 );
  xmaxg = 1.0 / sqrt ( realmin ( ) );
  xbnd = sqrt ( 50.0 );

  if ( x < 4.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R89SIFG - Fatal error!\n' );
    fprintf ( 1, '  Approximation invalid for X < 4.\n' );
    error ( 'R89SIFG - Fatal error!' )
  elseif ( x <= xbnd )
    value = csevl ( ( 1.0 ./ x ./ x - 0.04125 ) / 0.02125, f1cs, nf1 );
    f = ( 1.0 + value ) ./ x;
    value = csevl ( ( 1.0 ./ x ./ x - 0.04125 ) / 0.02125, g1cs, ng1 );
    g = ( 1.0 + value ) ./ x ./ x;
  elseif ( x <= xbig )
    value = csevl ( 100.0 ./ x ./ x - 1.0, f2cs, nf2 );
    f = ( 1.0 + value ) ./ x;
    value = csevl ( 100.0 ./ x ./ x - 1.0, g2cs, ng2 );
    g = ( 1.0 + value ) ./ x ./ x;
  else

    if ( x < xmaxf )
      f = 1.0 ./ x;
    else
      f = 0.0;
    end

    if ( x < xmaxg )
      g = 1.0 ./ x ./ x;
    else
      g = 0.0;
    end

  end

  return
end
