## ratint.m -  Octave file
##    begin                : Thu Oct 28 2003
##    copyright            : (C) 2003 by Marek Rychlik
##    email                : rychlik@u.arizona.edu
############################################################################
##                                                                         #
##   This program is free software; you can redistribute it and/or modify  #
##   it under the terms of the GNU General Public License as published by  #
##   the Free Software Foundation; either version 2 of the License, or     #
##   (at your option) any later version.                                   #
##                                                                         #
###########################################################################*

##
## SYNOPSIS: ratint(B, A, lo, hi, TOL)
##
## DESCRIPTION: Integrates any ratio of two polynomials B/A over the
## interval [lo, hi]. Returns the value of the integral.
##
## COMMENTS: The value is based on the partial fraction decomposition,
## and thus it is very accurate;it should exceed the accuracy of
## numerical quadrature.
##

function retval = ratint(B, A, lo, hi, TOL)
  ## Handle arguments
  if(nargin < 4)
    error("ratint(B, A, LO, HI, TOL), where TOL = 10^-3 if not \
	supplied.");
  endif
  if(nargin < 5)
    TOL = 1e-3;
  endif
  ## Find partial fractions
  [R, P, K, E] = residue(B, A, TOL);
  empty_list_elements_ok = 1;	# So that the next line does not pesture
				# us with a warning
  y = [];
  for x = [lo, hi];
    val = 0;
    ## Evaluate the following except when constant denominator
    if(length(E) > 0) 		
      E1 = (E ~= 1);		# These will be simple fractions
      E2 = (E == 1);		# These will be the logarithms
      ## Evaluate the simple fraction part
      if(any(E1))
	F1 = E(E1);
	R1 = R(E1);
	P1 = P(E1);
	## Prepare powers for integral
	--F1;
	val = val - sum((R1 ./ F1) ./ ((x - P1) .^ F1));
      endif
      ## Evaluate the logarithm part
      if(any(E2)) 
	## F2 = E(E2);
	R2 = R(E2);
	P2 = P(E2);
	val = val + sum(R2 .* log(x - P2));
      endif
    endif
    ## Evaluate the polynomial part
    if(length(K) > 0)
      val = val + polyval(polyinteg(K), x);
    endif
    ## Accumulate value
    y = [y, val];
  endfor
  ## Apply Fundamental Theorem of Calculus
  retval = y(2) - y(1);
endfunction

## ;;; Local Variables: ***
## ;;; mode: octave ***
## ;;; End: ***
