function value = ignpoi ( mu )

%*****************************************************************************80
%
%% IGNPOI generates a Poisson random deviate.
%
%  Discussion:
%
%    This procedure generates a single random deviate from a Poisson
%    distribution with given mean.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2014
%
%  Author:
%
%    Original FORTRAN77 version by Barry Brown, James Lovato.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Joachim Ahrens, Ulrich Dieter,
%    Computer Generation of Poisson Deviates
%    From Modified Normal Distributions,
%    ACM Transactions on Mathematical Software,
%    Volume 8, Number 2, June 1982, pages 163-179.
%
%  Parameters:
%
%    Input, real MU, the mean of the Poisson distribution
%    from which a random deviate is to be generated.
%
%    Output, integer VALUE, a random deviate from
%    the distribution.
%
  a0 = -0.5;
  a1 =  0.3333333;
  a2 = -0.2500068;
  a3 =  0.2000118;
  a4 = -0.1661269;
  a5 =  0.1421878;
  a6 = -0.1384794;
  a7 =  0.1250060;

  fact = [ 1.0, 1.0, 2.0, 6.0, 24.0, ...
    120.0, 720.0, 5040.0, 40320.0, 362880.0 ];
%
%  Start new table and calculate P0.
%
  if ( mu < 10.0 )

    m = max ( 1, floor ( mu ) );
    l = 0;
    p = exp ( - mu );
    q = p;
    p0 = p;
%
%  Uniform sample for inversion method.
%
    while ( true )

      u = r4_uni_01 ( );
      value = 0;

      if ( u <= p0 )
        return;
      end
%
%  Creation of new Poisson probabilities.
%
      for k = 1 : 35
        p = p * mu / k;
        q = q + p;
        if ( u <= q )
          value = k;
          return
        end
      end

    end

  else

    s = sqrt ( mu );
    d = 6.0 * mu * mu;
    l = floor ( mu - 1.1484 );
%
%  Normal sample.
%
    g = mu + s * snorm ( );

    if ( 0.0 <= g )

      value = floor ( g );
%
%  Immediate acceptance if large enough.
%
      if ( l <= value )
        return
      end
%
%  Squeeze acceptance.
%
      fk = value;
      difmuk = mu - fk;
      u = r4_uni_01 ( );

      if ( difmuk * difmuk * difmuk <= d * u )
        return
      end

    end
%
%  Preparation for steps P and Q.
%
    omega = 0.3989423 / s;
    b1 = 0.04166667 / mu;
    b2 = 0.3 * b1 * b1;
    c3 = 0.1428571 * b1 * b2;
    c2 = b2 - 15.0 * c3;
    c1 = b1 - 6.0 * b2 + 45.0 * c3;
    c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
    c = 0.1069 / mu;

    if ( 0.0 <= g )

      kflag = 0;

      if ( value < 10 )

        px = -mu;
        py = mu ^ value / fact(value+1);

      else

        del = 0.8333333E-01 / fk;
        del = del - 4.8 * del * del * del;
        v = difmuk / fk;

        if ( 0.25 < abs ( v ) )
          px = fk * log ( 1.0 + v ) - difmuk - del;
        else
          px = fk * v * v * ((((((( a7 ...
            * v + a6 ) ...
            * v + a5 ) ...
            * v + a4 ) ...
            * v + a3 ) ...
            * v + a2 ) ...
            * v + a1 ) ...
            * v + a0 ) - del;
        end

        py = 0.3989423 / sqrt ( fk );

      end

      x = ( 0.5 - difmuk ) / s;
      xx = x * x;
      fx = -0.5 * xx;
      fy = omega * ((( c3 * xx + c2 ) * xx + c1 ) * xx + c0 );

      if ( fy - u * fy <= py * exp ( px - fx ) )
        return
      end

    end
%
%  Exponential sample.
%
    while ( true )

      e = sexpo ( );
      u = 2.0 * r4_uni_01 ( ) - 1.0;
      if ( u < 0.0 )
        t = 1.8 - abs ( e );
      else
        t = 1.8 + abs ( e );
      end

      if ( t <= -0.6744 )
        continue
      end

      value = floor ( mu + s * t );
      fk = value;
      difmuk = mu - fk;

      kflag = 1;
%
%  Calculation of PX, PY, FX, FY.
%
      if ( value < 10 )

        px = -mu;
        py = mu ^ value / fact(value+1);

      else

        del = 0.8333333E-01 / fk;
        del = del - 4.8 * del * del * del;
        v = difmuk / fk;

        if ( 0.25 < abs ( v ) )
          px = fk * log ( 1.0 + v ) - difmuk - del;
        else
          px = fk * v * v * ((((((( a7 ...
            * v + a6 ) ...
            * v + a5 ) ...
            * v + a4 ) ...
            * v + a3 ) ...
            * v + a2 ) ...
            * v + a1 ) ...
            * v + a0 ) - del;
        end

        py = 0.3989423 / sqrt ( fk );

      end

      x = ( 0.5 - difmuk ) / s;
      xx = x * x;
      fx = -0.5 * xx;
      fy = omega * ((( c3 * xx + c2 ) * xx + c1 ) * xx + c0 );

      if ( kflag <= 0 )

        if ( fy - u * fy <= py * exp ( px - fx ) )
          return
        end

      else

        if ( c * abs ( u ) <= py * exp ( px + e ) - fy * exp ( fx + e ) )
          return
        end

      end

    end

  end

end
