function x = von_mises_sample ( a, b )

%*****************************************************************************80
%
%% von_mises_sample() samples the von Mises PDF.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    D J Best and N I Fisher,
%    Efficient Simulation of the von Mises Distribution,
%    Applied Statistics,
%    Volume 28, Number 2, pages 152-157.
%
%  Input:
%
%    real A, B, the parameters of the PDF.
%    -PI <= A <= PI,
%    0.0 < B.
%
%  Output:
%
%    real X, a sample of the PDF.
%
  tau = 1.0 + sqrt ( 1.0 + 4.0 * b * b );
  rho = ( tau - sqrt ( 2.0 * tau ) ) / ( 2.0 * b );
  r = ( 1.0 + rho * rho ) / ( 2.0 * rho );

  while ( true )

    u1 = rand ( 1, 1 );
    z = cos ( pi * u1 );
    f = ( 1.0 + r * z ) / ( r + z );
    c = b * ( r - f );

    u2 = rand ( 1, 1 );

    if ( u2 < c * ( 2.0 - c ) )
      break
    end

    if ( c <= log ( c / u2 ) + 1.0 )
      break
    end

  end

  u3 = rand ( 1, 1 );

  if ( u3 < 0.5 )
    x = a - acos ( f );
  else
    x = a + acos ( f );
  end

  return
end
