function P = lyngby_ks_prob(M, tol)

% lyngby_ks_prob       - Kolmogorov-Smirnov probability
%
%       function P = lyngby_ks_prob(M, tol)
%
%       Input: M     The Kolmogorov-Smirnov test value
%              tol   (optional) {eps} precision required for 
%                    the convergence. The default value is MATLAB's
%                    eps. 
%    
%       This function computes the function that enters 
%       into the calculation of the significance level in a 
%       Kolmogorov-Smirnov test.
%
%                      +oo
%                     .-- [    k-1         2  2  ]
%          Q  (x) = 2  >  [ (-1)   exp(-2.x .k ) ]
%           KS        '-- [                      ]
%                      k=1
%
%       The probability is calculated from an infinite series. If a
%       term in the series is lower than 'tol' the iteration will stop.
%       The algorithm will maximal use 100 terms.
%      
%       The probability is close to one the algorithm will use many
%       terms for the convergence. The user might want to lower the
%       value of 'tol' to get a faster (though less precise)
%       computation in the case where many data point have a
%       probability close to one.
%
%       Ref.: Numerical Recipes, 2nd ed. pp.623-628
%
%       See also: lyngby_ks_test, lyngby_ks_main

% cvs : $Id: lyngby_ks_prob.m,v 1.5 1997/09/04 13:47:27 fnielsen Exp $
%       $Revision: 1.5 $

    if (nargin < 1 | nargin > 2)
      error('lyngby_ks_test requires one or two arguments.') ;
    end

    if (nargin == 1)
      tol = eps ;
    else
      if (tol <= 0)
	error('''tol'' must be strictly positive.') ;
      end
    end

    [nl, nc] = size(M);
    if ((nl * nc) <= 0)
      error('First argument is empty or has wrong dimensions.');
    end

    n = nl * nc;
    X = reshape(M, n, 1);
    P = ones(n, 1);

    for i = find(X ~= 0)'
      coeff = 2;
      alpha = -2 * X(i, 1) * X(i, 1);
      sss   = 0;
      k     = 1;
      ui    = coeff * exp(alpha);
      while (( abs(ui) > tol * sss ) & ( k < 100 ))
	sss = sss + ui;
	k = k + 1;
	coeff = - coeff;
	ui = coeff * exp(alpha * k * k);
      end
      if (k == 100)
	% No converge
	P(i,1) = 1;
      else
	% It converges
	P(i,1) = round((sss + ui) / tol) * tol;
      end
    end

    P = reshape(P, nl, nc);













