function value = r4_dawson ( x )

%*****************************************************************************80
%
%% R4_DAWSON evaluates Dawson's integral of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the value of Dawson's integral at X.
%
  persistent daw2cs
  persistent dawacs
  persistent dawcs
  persistent ntdaw
  persistent ntdaw2
  persistent ntdawa
  persistent xbig
  persistent xmax
  persistent xsml

  if ( isempty ( ntdaw ) )

    daw2cs = [ ...
     -0.056886544105215527E+00, ...
     -0.31811346996168131E+00, ...
      0.20873845413642237E+00, ...
     -0.12475409913779131E+00, ...
      0.067869305186676777E+00, ...
     -0.033659144895270940E+00, ...
      0.015260781271987972E+00, ...
     -0.006348370962596214E+00, ...
      0.002432674092074852E+00, ...
     -0.000862195414910650E+00, ...
      0.000283765733363216E+00, ...
     -0.000087057549874170E+00, ...
      0.000024986849985481E+00, ...
     -0.000006731928676416E+00, ...
      0.000001707857878557E+00, ...
     -0.000000409175512264E+00, ...
      0.000000092828292216E+00, ...
     -0.000000019991403610E+00, ...
      0.000000004096349064E+00, ...
     -0.000000000800324095E+00, ...
      0.000000000149385031E+00, ...
     -0.000000000026687999E+00, ...
      0.000000000004571221E+00, ...
     -0.000000000000751873E+00, ...
      0.000000000000118931E+00, ...
     -0.000000000000018116E+00, ...
      0.000000000000002661E+00, ...
     -0.000000000000000377E+00, ...
      0.000000000000000051E+00 ]';
    dawacs = [ ...
      0.01690485637765704E+00, ...
      0.00868325227840695E+00, ...
      0.00024248640424177E+00, ...
      0.00001261182399572E+00, ...
      0.00000106645331463E+00, ...
      0.00000013581597947E+00, ...
      0.00000002171042356E+00, ...
      0.00000000286701050E+00, ...
     -0.00000000019013363E+00, ...
     -0.00000000030977804E+00, ...
     -0.00000000010294148E+00, ...
     -0.00000000000626035E+00, ...
      0.00000000000856313E+00, ...
      0.00000000000303304E+00, ...
     -0.00000000000025236E+00, ...
     -0.00000000000042106E+00, ...
     -0.00000000000004431E+00, ...
      0.00000000000004911E+00, ...
      0.00000000000001235E+00, ...
     -0.00000000000000578E+00, ...
     -0.00000000000000228E+00, ...
      0.00000000000000076E+00, ...
      0.00000000000000038E+00, ...
     -0.00000000000000011E+00, ...
     -0.00000000000000006E+00, ...
      0.00000000000000002E+00 ]';
    dawcs = [ ...
     -0.006351734375145949E+00, ...
     -0.22940714796773869E+00, ...
      0.022130500939084764E+00, ...
     -0.001549265453892985E+00, ...
      0.000084973277156849E+00, ...
     -0.000003828266270972E+00, ...
      0.000000146285480625E+00, ...
     -0.000000004851982381E+00, ...
      0.000000000142146357E+00, ...
     -0.000000000003728836E+00, ...
      0.000000000000088549E+00, ...
     -0.000000000000001920E+00, ...
      0.000000000000000038E+00 ]';

    eps = r4_mach ( 3 );
    ntdaw  = r4_inits ( dawcs,  13, 0.1 * eps );
    ntdaw2 = r4_inits ( daw2cs, 29, 0.1 * eps );
    ntdawa = r4_inits ( dawacs, 26, 0.1 * eps );

    xsml = sqrt ( 1.5 * eps );
    xbig = sqrt ( 0.5 / eps );
    xmax = exp ( min ( - log ( 2.0 * r4_mach ( 1 ) ), ...
      log ( r4_mach ( 2 ) ) ) - 0.01 );

  end

  y = abs ( x );

  if ( y <= xsml )
    value = x;
  elseif ( y <= 1.0 )
    value = x * ( 0.75 + r4_csevl ( 2.0 * y * y - 1.0, dawcs, ntdaw ) );
  elseif ( y <= 4.0 )
    value = x * ( 0.25 + r4_csevl ( 0.125 * y * y - 1.0, daw2cs, ntdaw2 ) );
  elseif ( y < xbig )
    value = ( 0.5 + r4_csevl ( 32.0 / y / y - 1.0, dawacs, ntdawa ) ) / x;
  elseif ( y <= xmax )
    value = 0.5 / x;
  else
    value = 0.0;
  end

  return
end