!********************************************************
!*    Collection of Chebyshev approximation routines    *
!* ---------------------------------------------------- *
!* REFERENCE: "Numerical Recipes, The Art of Scientific *
!*             Computing By W.H. Press, B.P. Flannery,  *
!*             S.A. Teukolsky and W.T. Vetterling,      *
!*             Cambridge University Press, 1986"        *
!*             [BIBLI 08].                              *             
!*                                                      *
!*                F90 Release 1.1 By J-P Moreau, Paris. *
!*                        (www.jpmoreau.fr)             *
!* ---------------------------------------------------- *
!* Release 1.1: added subroutines CHINT and CHDER.      *
!********************************************************

subroutine CHEBFT(A,B,C,N)
real*8 A,B,C(N)
external FUNC
!********************************************************
!* Chebyshev fit: Given a real function FUNC(X), lower  *
!* and upper limits of the interval [A,B] for X, and a  *
!* maximum degree N, this routine computes the N Cheby- *
!* shev coefficients Ck, such that FUNC(X) is approxima-*
!* ted by:  N                                           *
!*         [Sum Ck Tk-1(Y)] - C1/2, where X and Y are   *
!*         k=1                                          *
!* related by:     Y = (X - 1/2(A+B)) / (1/2(B-A))      *
!* This routine is to be used with moderately large N   *
!* (e.g. 30 or 50), the array of C's subsequently to be *
!* truncated at the smaller value m such that Cm+1 and  *
!* subsequent elements are negligible.                  *
!********************************************************
parameter(NMAX=50,HALF=0.5d0,TWO=2.d0,ZERO=0.d0)
real*8 PI,SUM,F(NMAX)
real*8 BMA,BPA,FAC, Y
real*8 FUNC
PI=4.d0*datan(1.d0)
BMA=HALF*(B-A); BPA=HALF*(B+A)
do K=1,N
  Y=DCOS(PI*(K-HALF)/N)
  F(K)=FUNC(Y*BMA+BPA)
end do
FAC=TWO/N
do J=1,N
  SUM=ZERO
  do K=1,N
    SUM=SUM+F(K)*DCOS((PI*(J-1))*((K-HALF)/N))
  end do
  C(J)=FAC*SUM
end do
return
end

real*8 function CHEBEV(A,B,C,M,X)
!**********************************************************
!* Chebyshev evaluation: All arguments are input. C is an *
!* array of Chebyshev coefficients, of length M, the first*
!* M elements of Coutput from subroutine CHEBFT (which    *
!* must have been called with the same A and B). The Che- *
!* byshev polynomial is evaluated at a point Y determined *
!* from X, A and B, and the result FUNC(X) is returned as *
!* the function value.                                    *
!**********************************************************
parameter(HALF=0.5d0,TWO=2.d0,ZERO=0.d0)
real*8 A,B,C(M),X
real*8 D,DD,SV,Y,Y2
if ((X-A)*(X-B).gt.ZERO) pause ' X not in range.'
D=ZERO; DD=ZERO
Y=(TWO*X-A-B)/(B-A)  !change of variable
Y2=TWO*Y
do J=M,2,-1
  SV=D
  D=Y2*D-DD+C(J)
  DD=SV
end do
CHEBEV=Y*D-DD+HALF*C(1)
return
end

subroutine CHINT(A,B,C,CINT,N)
!**********************************************************
!* Given A,B,C, as output from routine CHEBFT, and given  *
!* N, the desired degree of approximation (length of C to *
!* be used), this routine returns the array CINT, the Che-*
!* byshev coefficients of the integral of the function    *
!* whose coefficients are C. The constant of integration  *
!* is set so that the integral vanishes at A.             *
!**********************************************************
real*8 A,B,C(N),CINT(N)
parameter(ONE=1.d0,QUART=0.25d0,TWO=2.d0,ZERO=0.d0)
real*8 CON,FAC,SUM
CON=QUART*(B-A)
SUM=ZERO
FAC=ONE
do J=2, N-1
  CINT(J)=CON*(C(J-1)-C(J+1))/(J-1)
  SUM=SUM+FAC*CINT(J)
  FAC=-FAC
end do
CINT(N)=CON*C(N-1)/(N-1)
SUM=SUM+FAC*CINT(N)
CINT(1)=TWO*SUM  !set the constant of integration
return
end

subroutine CHDER(A,B,C,CDER,N)
!**********************************************************
!* Given A,B,C, as output from routine CHEBFT, and given  *
!* N, the desired degree of approximation (length of C to *
!* be used), this routine returns the array CDER, the Che-*
!* byshev coefficients of the derivative of the function  *
!* whose coefficients are C.                              *
!********************************************************** 
real*8 A,B,C(N),CDER(N)
parameter(TWO=2.d0,ZERO=0.d0)
real*8 CON
CDER(N)=ZERO
CDER(N-1)=TWO*(N-1)*C(N)
if (N.ge.3) then
  do J=N-2, 1, -1
    CDER(J)=CDER(J+2)+TWO*J*C(J+1)
  end do
end if
CON=TWO/(B-A)
do J=1,N                !normalize to interval B - A
  CDER(J)=CDER(J)*CON
end do
return
end

!end of file Chebyshe.f90