function x = lyngby_idf_f(P, r1, r2)

% lyngby_idf_f         - Inverse cumulative F distribution.
% 
%      function x =  lyngby_idf_f(P, r1, r2);
%
%      Input:  p    The requested probability P(T<p) where F
%                   is F-distributed with ...
%              r1   Degrees of freedom in the nominator.
%              r2   Degrees of freedom in the denominator.
%
%      Output: x  The corresponding value if F.
%
%      This function will return the inverse of the cumulative
%      Fisher's F probability distribution. 
%
%      Note that this function is NOT vectorized !
%
%      See also: lyngby_cdf_f, lyngby_pdf_f, lyngby_idf_t,
%                lyngby_idf_chi2

% cvs : $Id: lyngby_idf_f.m,v 1.2 1999/09/16 13:07:06 fnielsen Exp $ 
%       $Revision: 1.2 $    


   if nargin < 3, 
     error('This function requires three inputs.'); 
   end

   if prod(size(P)) ~= 1
     error('P should only contain one value');
   end

   evalType = 4;
   

   if evalType == 1
     [x, fval, exitflag] = ...
	 fzero(sprintf('lyngby_cdf_f(x,%i,%i)-%f',r1,r2,P), ...
	 4, optimset('Display', 'off'));
   
   elseif evalType == 2

     % This one does not work in Matlab 4
     x = fzero(sprintf('lyngby_cdf_f(x,%i,%i)-%f',r1,r2,P), ...
	 [1e-6 1000],eps, 1);
     
   elseif evalType == 3
     
     % This one does not work in Matlab 4, and usually not in Matlab 5
     x = fzero(sprintf('lyngby_cdf_f(x,%i,%i)-%f',r1,r2,P)', 3.4, ...
	 eps, 1); 

   else
     
     % Bisection
     x1 = 0.001;
     x2 = 1000;
     for e = 1:100
       x = (x1+x2)/2;
       Piter = lyngby_cdf_f(x, r1, r2);
       %     [x x1 x2 P Piter]
       % disp([x Piter]), keyboard;
       if Piter < P
	 x1 = x;
       else
	 x2 = x;
       end
     end
     
   end
















