function pdf = fisher_pdf ( x, kappa, mu )

%*****************************************************************************80
%
%% FISHER_PDF evaluates the Fisher PDF.
%
%  Discussion:
%
%    The formulat for the PDF is:
%
%      PDF(KAPPA,MU;X) = C(KAPPA) * exp ( KAPPA * MU' * X )
%
%    where:
%
%      0 <= KAPPA is the concentration parameter,
%      MU is a point on the unit sphere, the mean direction,
%      X is any point on the unit sphere,
%      and C(KAPPA) is a normalization factor:
%
%      C(KAPPA) = sqrt ( KAPPA ) / ( ( 2 * pi )**(3/2) * I(0.5,KAPPA) )
%
%    where
%
%      I(nu,X) is the Bessel function of order NU and argument X.
%
%    For a fixed value of MU, the value of KAPPA determines the
%    tendency of sample points to tend to be near MU.  In particular,
%    KAPPA = 0 corresponds to a uniform distribution of points on the
%    sphere, but as KAPPA increases, the sample points will tend to
%    cluster more closely to MU.
%
%    The Fisher distribution for points on the unit sphere is
%    analogous to the normal distribution of points on a line,
%    and, more precisely, to the von Mises distribution on a circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Kanti Mardia, Peter Jupp,
%    Directional Statistics,
%    Wiley, 2000,
%    LC: QA276.M335
%
%  Parameters:
%
%    Input, real X(3), the argument of the PDF.
%    X should have unit Euclidean norm, but this routine will
%    automatically work with a normalized version of X.
%
%    Input, real KAPPA, the concentration parameter.
%    0 <= KAPPA is required.
%
%    Input, real MU(3), the mean direction.
%    MU should have unit Euclidean norm, but this routine will
%    automatically work with a normalized version of MU.
%
%    Output, real PDF, the value of the PDF.
%

%
%  Force column-vector shape.
%
  x = x ( : );
  mu = mu ( : );

  nb = 1;

  if ( kappa < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FISHER_PDF - Fatal error!\n' );
    fprintf ( 1, '  KAPPA must be nonnegative.\n' );
    fprintf ( 1, '  Input KAPPA = %f\n', kappa );
    error ( 'FISHER_PDF - Fatal error!' );
  end

  if ( kappa == 0.0 )
    pdf = 1.0 / ( 4.0 * pi );
    return;
  end

  alpha = 0.5;
  ize = 1;

  [ b, ncalc ] = ribesl ( kappa, alpha, nb, ize );

  cf = sqrt ( kappa ) / ( sqrt ( ( 2.0 * pi )^3 ) * b(1) );

  mu_norm = norm ( mu );

  if ( mu_norm == 0.0 )
    pdf = cf;
    return
  end

  x_norm = norm ( x );

  if ( x_norm == 0.0 )
    pdf = cf;
    return
  end

  arg = kappa * ( x(1:3)' * mu(1:3) ) / ( x_norm * mu_norm );

  pdf = cf * exp ( arg );

  return
end