C+
C NAME:
C	POLYNOMIAL_EXP
C PURPOSE:
C	Given the coefficients of its expansion in some set of (orthogonal)
C	polynomials, evaluate the value of a polynomial in point X
C CATEGORY:
C	Math: orthogonal polynomials
C CALLING SEQUENCE:
	function POLYNOMIAL_EXP(A,B,X,A0,B0,N,PN,POLY)
C INPUTS:
C	A,B	real	(read-only); interval [A,B]; see PROCEDURE
C	X	real	(read-only); X-value in [A,B] where function is evaluated
C	A0,B0	real	(read-only); interval [A0,B0]; see PROCEDURE
C	N	integer	(read-only); # terms in polynomial expansion
C	PN(N)	real	(read-only); expansion coefficients
C	POLY	real	externally declared function, used to evaluate polynomials
C OUTPUTS:
C	POLYNOMIAL_EXP	real	value of polynomial in point X
C RESTRICTIONS:
C >	The number of terms in the expansion, N, must be larger then zero.
C >	A and B must be unequal.
C PROCEDURE:
C >	The function POLY must be declared external in the calling program
C	The call to POLY has the form F = POLY(XS,I); I is the degree of the
C	polynomial to be evaluated; XS = A0+(B0-A0)*(X-A)/(B-A) (i.e. XS is the
C	value corresponding to X after the interval [A,B] is rescaled to 
C	[A0,B0]).
C >	The input value of X is assumed to be in the same units as interval
C	[A,B]. The interval [A,B] is mapped to [A0,B0]; X is correspondingly
C	rescaled to XS=A0+(B0-A0)*(X-A)/(B-A). The expansion is evaluated in XS.
C >	The polynomial expansion is
C	VALUE(X) = SUM(i=1,N) { PN(I)*POLY_i-1(XS) }, i.e. PN(I) is the
C	coefficient of the polynomial of degree I-1
C >	The polynomials POLY are evaluated in function POLY (see EXTERNAL)
C MODIFICATION HISTORY:
C	DEC-1991; Paul Hick (UCSD)
C-
	    real	A
	    real	B
	    real	X
	    real	A0
	    real	B0
	    integer	N
	    real	PN(N)
	    real	POLY

	POLYNOMIAL_EXP = BadR4()

	if (N .gt. 0 .and. B .ne. A) then
	    XS = A0+(B0-A0)/(B-A)*(X-A)
	    POLYNOMIAL_EXP = 0.
	    do I=0,N-1
		POLYNOMIAL_EXP = POLYNOMIAL_EXP+PN(I+1)*POLY(XS,I) ! Expansion
	    end do
	end if

	return
	end