!	*********************************************************************************
!	*																				*
!	* This FORTRAN90/95 module is used to perform Chebyshev regression of			*
!	* functions up to four dimensions. The file chebydoc.tex (or chebydoc.pdf)		*
!	* contains a detailed description of the subroutines and the algorithm used.	*
!	* Note: The module can be modified to allow for approximation up to	seven		*
!	*	    dimensions.																*
!	*																				*
!	* Written by Constantino Hevia, June 2006.										*
!	*																				* 
!	*********************************************************************************

MODULE chebyshev_regression
IMPLICIT NONE

! Set precision
INTEGER, PARAMETER, PRIVATE :: wp = SELECTED_REAL_KIND(15)

CONTAINS

	SUBROUTINE chebnodes( m, x, lb, ub )

	! PURPOSE
	! To compute the chebyshev nodes on an interval [lb,ub]. The input variables
	! lb and ub are optional. If they are not present, then the standard nodes on [-1,1]
	! are computed.
	!
	! INPUTS
	! m		: number of knots
	! lb	: lower bound of approximation interval ( OPTIONAL ARGUMENT )
	! ub	: upper bound of approximation interval ( OPTIONAL ARGUMENT )
	!
	! OUTPUT
	! x		: array of dimension m with the chebyshev nodes
	! 
	! ROUTINES USED
	! none

	IMPLICIT NONE
	INTEGER, INTENT(in) :: m
	REAL(wp), INTENT(in), OPTIONAL :: lb, ub
	REAL(wp), INTENT(out) :: x( m )

	! Local variables
	INTEGER ::i				! loop index
	REAL(wp), PARAMETER		:: pi = 3.141592653589793238462643383279502884197_wp	! number pi

	! Compute the 'm' Chebyshev's nodes in [-1,1]
	DO i = 1, m
		x( i ) = - COS( pi*( 2.0_wp*i - 1.0_wp ) / ( 2.0_wp*DBLE(m) ) )
	END DO

	! Check if need to move nodes to [lb,ub]
	IF ( PRESENT( lb ) .and. PRESENT( ub ) ) THEN
		x = ( x + 1.0_wp ) * 0.5_wp * ( ub - lb ) + lb
	END IF
	
	END SUBROUTINE chebnodes


	SUBROUTINE chebbas( m, n , basis, z )

	! 	PURPOSE:
	!	To compute the Chebyshev basis matrix in one dimension
	!
	!	INPUTS:
	!	n		: degree of chebyshev polynomial (i.e. higher order term is x ** n )
	!	m		: number of nodes ( m >= n + 1 )
	!	
	!	OUTPUTS:
	!	basis	: Chebyshev's m by n + 1 basis matrix in one dimension	
	!	z		: nodes in [-1,1]
	!

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: m						! Number of knots
	INTEGER, INTENT(in)		:: n						! Degree of the polynomial
	REAL(wp), INTENT(out)	:: basis( m, 0 : n )		! Chebyshev's basis matrix
	REAL(wp), INTENT(out)	:: z( m )					! Chebyshev nodes in [-1,1]
			
	!Local variables
	INTEGER ::i

	CALL chebnodes( m, z )	! Compute the Chebyshev nodes in [-1,1]

	basis( :, 0 ) = 1.0_wp	! First column of the basis matrix
	basis( :, 1 ) = z		! Second column of the basis matrix
	DO i = 2, n  			! Additional columns of the basis matrix ( only done if n>=2 )
		basis( :, i ) = 2.0_wp * z * basis( : , i - 1 ) - basis( : , i - 2 )
	END DO
	
	END SUBROUTINE chebbas


	SUBROUTINE chebco1d( y, m, n, coeff )

	! PURPOSE:
	! To compute the Chebyshev's coefficients of a polynomial approximation of one dimension.
	! Given a function f( x1 ) find the Chebyshev's approximation of degree n.
	! This follows algorithm 6.2 in Judd.
	!
	! INPUTS:
	! m		: integer with the number of evaluation nodes
	! n		: integer with the degree of the approximation ( m>=n+1 )
	! y		: array of size m with the function's values at each node
	!
	! OUTPUT:
	! coeff = array of size (0:n) with the Chebyshev's coefficients

	IMPLICIT NONE
	INTEGER, INTENT(in) :: m, n
	REAL(wp), INTENT(in) :: y( m )
	REAL(wp), INTENT(out) :: coeff( 0:n )

	! Local variables
	REAL(wp) :: z( m ), t( m, 0:n ), const
	INTEGER ::i

	! check dimensions
	IF ( m <= n ) THEN
		WRITE(*,*) ' Error in subroutine chebco1d in chebyshev regression module.'
		WRITE(*,*) ' m <= n. '
		STOP
	END IF
	
	! Compute the basis and knots
	CALL chebbas( m, n , t, z )
	
	! Compute the Chebyshev coefficients
	coeff( 0 ) = SUM( y ) / DBLE( m )
	const = 2.0_wp / DBLE( m )
	DO i = 1, n
		coeff( i ) = const * DOT_PRODUCT( y, t( :, i ) ) 
	END DO

	END SUBROUTINE chebco1d



	SUBROUTINE chebco2d( y, m1, m2, n1, n2, coeff )

	! PURPOSE:
	! To compute the Chebyshev's coefficients of a polynomial approximation of two dimensions.
	! Given a function f( x1, x2 ) find the Chebyshev's approximation of degree n1 * n2
	! This subroutine follows algorithm 6.4 in Judd.
	!
	! INPUTS
	! m1	: number of nodes in dimension 1
	! m2	: number of nodes in dimension 2
	! n1	: degree of polynomial in dimension 1 ( m1 >= n1+1 )
	! n2	: degree of polynomial in dimension 2 ( m2 >= n2+1 )
	! y		: array of size ( m1, m2 ) with the function's values at each node
	!
	! OUTPUT
	! coeff : array of size( 0:n1, 0:n2 ) with the Chebyshev's coefficients

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: m1, m2, n1, n2
	REAL(wp), INTENT(in)	:: y( m1, m2 )
	REAL(wp), INTENT(out)	:: coeff( 0:n1, 0:n2 )

	! Local variables
	REAL(wp) :: z1( m1 ),z2( m2 ),t1( m1, 0:n1 ),t2( m2, 0:n2 ),ct1( 0:n1 ),ct2( 0:n2 ),temp
	INTEGER  :: i1, i2, j1, j2

	! Check degree and number of nodes
	IF ( (m1 <= n1) .or. (m2 <= n2) ) THEN
		WRITE( *, * ) ' Error in subroutine chebco2d in chebyshev regression module. '
		WRITE( *, * ) ' Either m1<=n1 or m2<=n2. '
		STOP
	END IF

	! Compute the basis and regression nodes in each dimension
	CALL chebbas( m1, n1 , t1, z1 )		! Dimension 1
	CALL chebbas( m2, n2 , t2, z2 )		! Dimension 2

	! Define constants - these come from the discrete orthogonality of chebyshev polynomials
	ct1( 0 )	 = 1.0_wp / DBLE( m1 )
	ct1( 1: n1 ) = 2.0_wp / DBLE( m1 )
	ct2( 0 )     = 1.0_wp / DBLE( m2 )
	ct2( 1: n2 ) = 2.0_wp / DBLE( m2 )

	! Compute Chebyshev coefficients
	DO i1 = 0, n1
		DO i2 = 0, n2
			temp = 0.0_wp
			DO j1 = 1, m1
				DO j2 = 1, m2
					temp = temp + y( j1, j2 ) * t1( j1, i1 ) * t2( j2, i2 )
				END DO
			END DO
			coeff( i1, i2 ) = ct1( i1 ) * ct2( i2 ) * temp
		END DO
	END DO

	END SUBROUTINE chebco2d

	

	SUBROUTINE chebco3d( y, m1, m2, m3, n1, n2, n3, coeff )

	! PURPOSE
	! To compute the Chebyshev's coefficients of a polynomial approximation of three dimensions.
	! Given a function f( x1, x2, x3 ) find the Chebyshev's approximation of order n1*n2*n3
	! This subroutine follows algorithm 6.4 in Judd.
	!
	! INPUTS
	! m1	: number of nodes in dimension 1
	! m2	: number of nodes in dimension 2
	! m3	: number of nodes in dimension 3
	! n1	: degree of polynomial in dimension 1 ( m1 >= n1+1 ) 
	! n2	: degree of polynomial in dimension 2 ( m2 >= n2+1 )
	! n3	: degree of polynomial in dimension 3 ( m3 >= n3+1 )
	! y		: array of size( m1, m2, m3 ) with the function's value at each node.
	!
	! OUTPUT:
	! coeff	: array of size( 0:n1, 0:n2, 0:n3 ) with the Chebyshev's coefficients
	
	IMPLICIT NONE
	INTEGER, INTENT(in)		:: m1, m2, m3, n1, n2, n3
	REAL(wp), INTENT(in)	:: y( m1, m2, m3 )
	REAL(wp), INTENT(out)	:: coeff( 0:n1, 0:n2, 0:n3 )

	! Local variables
	REAL(wp) :: z1( m1 ), z2( m2 ), z3( m3 ), t1( m1, 0:n1 ), t2( m2, 0:n2 ), t3( m3, 0:n3),	&
			    & temp, ct1( 0:n1 ), ct2( 0:n2 ), ct3( 0:n3 )
	INTEGER  :: i1, i2, i3, j1, j2, j3

	! Check degree and number of nodes
	IF ( (m1 <= n1) .or. (m2 <= n2) .or. (m3<=n3) ) THEN
		WRITE( *, * ) ' Error in subroutine chebco3d in chebyshev regression module. '
		WRITE( *, * ) ' Either m1<=n1, m2<=n2 or m3<=n3. '
		STOP
	END IF

	! Compute the basis and nodes in each dimension
	CALL chebbas( m1, n1, t1, z1 )		! Dimension 1
	CALL chebbas( m2, n2, t2, z2 )		! Dimension 2
	CALL chebbas( m3, n3, t3, z3 )		! Dimension 3

	! Define constants - these come from the discrete orthogonality of chebyshev polynomials
	ct1( 0 )	 = 1.0_wp / DBLE( m1 )
	ct1( 1: n1 ) = 2.0_wp / DBLE( m1 )
	ct2( 0 )     = 1.0_wp / DBLE( m2 )
	ct2( 1: n2 ) = 2.0_wp / DBLE( m2 )
	ct3( 0 )     = 1.0_wp / DBLE( m3 )
	ct3( 1: n3 ) = 2.0_wp / DBLE( m3 )

	! Compute Chebyshev's coefficients
	DO i1 = 0, n1
		DO i2 = 0, n2
			DO i3 = 0, n3
				temp = 0.0_wp
				DO j1 = 1, m1
					DO j2 = 1, m2
						DO j3 = 1, m3
							temp = temp + y(j1,j2,j3)*t1(j1,i1)*t2(j2,i2)*t3(j3,i3)
						END DO
					END DO
				END DO
				coeff(i1,i2,i3) = ct1( i1 ) * ct2( i2 ) * ct3( i3 ) * temp
			END DO
		END DO
	END DO

	END SUBROUTINE chebco3d

	
	SUBROUTINE chebco4d( y, m1, m2, m3, m4, n1, n2, n3, n4, coeff )

	! PURPOSE
	! To compute the Chebyshev's coefficients of a polynomial approximation of four dimensions.
	! Given a function f( x1, x2, x3, x4 ) find the Chebyshev's approximation of order n1*n2*n3*n4
	! This subroutine follows algorithm 6.4 in Judd.
	! Warning: the tensor product approach can make this routine very time intensive for high
	!		   degrees of interpolation 
	!
	! INPUTS
	! m1	: number of nodes in dimension 1
	! m2	: number of nodes in dimension 2
	! m3	: number of nodes in dimension 3
	! m4	: number of nodes in dimension 4
	! n1	: degree of polynomial in dimension 1 ( m1 >= n1+1 )
	! n2	: degree of polynomial in dimension 2 ( m2 >= n2+1 )
	! n3	: degree of polynomial in dimension 3 ( m3 >= n3+1 )
	! n4	: degree of polynomial in dimension 4 ( m4 >= n4+1 )
	! y		: array of size( m1, m2, m3, m4 ) with the function's value at each node.
	!
	! OUTPUT:
	! coeff	: array of size( 0:n1, 0:n2, 0:n3, 0:n4 ) with the Chebyshev's coefficients

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: m1, m2, m3, m4, n1, n2, n3, n4
	REAL(wp), INTENT(in)	:: y( m1, m2, m3, m4 )
	REAL(wp), INTENT(out)	:: coeff( 0:n1, 0:n2, 0:n3, 0:n4 )


	! Local variables
	REAL(wp) :: z1( m1 ), z2( m2 ), z3( m3 ), z4( m4 ), t1( m1, 0:n1 ), t2( m2, 0:n2 ),			&
				& t3( m3, 0:n3 ), t4( m4, 0:n4 ), temp, ct1( 0:n1 ), ct2( 0:n2 ), ct3( 0:n3 ),	&
				& ct4( 0:n4 )
	INTEGER  :: i1, i2, i3, i4, j1, j2, j3, j4

	! Check degree and number of nodes
	IF ( (m1 <= n1) .or. (m2 <= n2) .or. (m3<=n3) .or. (m4<=n4) ) THEN
		WRITE( *, * ) ' Error in subroutine chebco3d in chebyshev regression module. '
		WRITE( *, * ) ' Either m1<=n1, m2<=n2, m3<=n3 or m4<=n4. '
		STOP
	END IF


	! Compute the basis and nodes in each dimension
	CALL chebbas( m1, n1, t1, z1 )		! Dimension 1
	CALL chebbas( m2, n2, t2, z2 )		! Dimension 2
	CALL chebbas( m3, n3, t3, z3 )		! Dimension 3
	CALL chebbas( m4, n4, t4, z4 )		! Dimension 4

	! Define constants - these come from the discrete orthogonality of chebyshev polynomials
	ct1( 0 )	 = 1.0_wp / DBLE( m1 )
	ct1( 1: n1 ) = 2.0_wp / DBLE( m1 )
	ct2( 0 )     = 1.0_wp / DBLE( m2 )
	ct2( 1: n2 ) = 2.0_wp / DBLE( m2 )
	ct3( 0 )     = 1.0_wp / DBLE( m3 )
	ct3( 1: n3 ) = 2.0_wp / DBLE( m3 )
	ct4( 0 )     = 1.0_wp / DBLE( m4 )
	ct4( 1: n4 ) = 2.0_wp / DBLE( m4 )


	! Compute Chebyshev's coefficients
	DO i1 = 0, n1
		DO i2 = 0, n2
			DO i3 = 0, n3
				DO i4 = 0, n4
					temp = 0.0_wp
					DO j1 = 1, m1
						DO j2 = 1, m2
							DO j3 = 1, m3
								DO j4 = 1, m4
									temp = temp + y(j1,j2,j3,j4)*t1(j1,i1)*t2(j2,i2)*t3(j3,i3)*	&
									&	   t4(j4,i4)
								END DO
							END DO
						END DO
					END DO
					coeff( i1, i2, i3, i4 ) = ct1( i1 )*ct2( i2 )*ct3( i3 )*ct4( i4 )*temp
				END DO
			END DO
		END DO
	END DO

	END SUBROUTINE chebco4d


	! THE FOLLOWING FUNCTIONS EVALUATE THE CHEBYSHEV'S POLYNOMIALS:


	FUNCTION chebev1d( n, lb, ub, x, coeff )

	! PURPOSE
	! To evaluate the 1-dimensional chebyshev approximate function at the value x. This function
	! uses Clenshaw algorithm, which is an efficient way of evaluating 1 dimensional Chebyshev 
	! polynomials. ( Mason and Handscomb, "Chebyshev Polynomials", 2003, section 2.4. )
	!
	! INPUTS
	! n		: degree of polynomial.
	! lb	: lower bound.
	! ub	: upper bound.
	! x		: evaluation point that must belong to [lb,ub].
	! coeff : array of size (0:n) witht the Chebyshev's coefficients
	!
	! OUTPUT
	! chebev1d	: value of the approximate function at x

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n
	REAL(wp), INTENT(in)	:: x, coeff( 0:n ), lb, ub
	REAL(wp) :: chebev1d

	! Local variables
	REAL(wp) :: z, b( 0: n+2 )
	INTEGER  :: i

	! Transform x to the interval [-1,1]
	z = 2.0_wp * ( x - lb )/( ub - lb ) - 1.0_wp

	! Initialize vector b
	b = 0.0_wp

	! Compute Clenshaw recursion
	DO i = n, 0, -1
		b( i ) = 2.0_wp * z * b( i + 1 ) - b( i + 2 ) + coeff( i )
	END DO

	! Evaluate chebyshev polynomial
	chebev1d = b( 0 ) - z * b( 1 )

	END FUNCTION chebev1d


	FUNCTION chebev2d( n1, n2, lb1, ub1, lb2, ub2, x1, x2, coeff )
	! PURPOSE
	! To evaluate the 2-dimensional chebyshev approximate function at the values (x1,x2). This 
	! function does not use Clenshaw recursion. In tests, this is more efficient than chebev2d2
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2) with the chebyshev's coefficients.
	!
	! OUTPUT
	! chebev2d	: value of the approximate function at (x1,x2)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2
	REAL(wp), INTENT(in)	:: x1, x2, coeff( 0:n1, 0:n2 ), lb1, ub1, lb2, ub2
	REAL(wp) :: chebev2d

	! Local variables
	REAL(wp) :: z, t1( 0:n1 ), t2( 0:n2 )
	INTEGER  :: i, j

	! 1) Evaluate the Chebyshev polynomials at the points x1 and x2
	! First dimension
	z = 2.0_wp*( x1 - lb1 ) /( ub1 - lb1 ) - 1.0_wp
	t1( 0 ) = 1.0_wp
	t1( 1 ) = z
	DO i = 2, n1
		t1( i ) = 2.0_wp * z * t1( i - 1 ) - t1( i - 2 )
	END DO

	! Second dimension
	z = 2.0_wp*( x2 - lb2 ) /( ub2 - lb2 ) - 1.0_wp
	t2( 0 ) = 1.0_wp
	t2( 1 ) = z
	DO i = 2, n2
		t2( i ) = 2.0_wp * z * t2( i - 1 ) - t2( i - 2 )
	END DO

	! 2) Evaluate the approximation
	chebev2d = 0.0_wp
	DO i = 0, n1
		DO j = 0, n2
			chebev2d = chebev2d + coeff( i, j ) * t1( i ) * t2( j )
		END DO
	END DO

	END FUNCTION chebev2d


	FUNCTION chebev2d2( n1, n2, lb1, ub1, lb2, ub2, x1, x2, coeff )

	! PURPOSE
	! To evaluate the 2-dimensional chebyshev approximate function at the values (x1,x2). This 
	! function uses Clenshaw algorithm. The one dimensional algorithm is extended to two 
	! dimensions (consult my notes). In tests, this is less efficient than chebev2d
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2) with the chebyshev's coefficients.
	!
	! OUTPUT
	! chebev2d2	: value of the approximate function at (x1,x2)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2
	REAL(wp), INTENT(in)	:: x1, x2, coeff( 0:n1, 0:n2 ), lb1, ub1, lb2, ub2
	REAL(wp) :: chebev2d2

	! Local variables
	REAL(wp) :: z, a( 0:n1 ), b( 0 : n2 + 2 )
	INTEGER  :: i,j

	! Transform x2 to the interval [-1,1]
	z = 2.0_wp * ( x2 - lb2 )/( ub2 - lb2 ) - 1.0_wp

	! Reduce evaluation to a 1 dimensional problem by summing along second dimension, using
	! Clenshaw recursions along second dimension.
	DO i = 0, n1
		b = 0.0_wp
		DO j = n2, 0, -1		
			b( j ) = 2.0_wp * z * b( j + 1 ) - b( j + 2 ) + coeff( i, j )
		END DO
		a( i ) = b( 0 ) - z * b( 1 )
	END DO

	! Evaluate Chebyshev aproximation as a 1 dimensional problem
	chebev2d2 = chebev1d( n1, lb1, ub1, x1, a )

	END FUNCTION chebev2d2


	FUNCTION chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3, coeff )

	! PURPOSE
	! To evaluate the 3-dimensional chebyshev approximate function at the values (x1,x2,x3)
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3) with the chebyshev's coefficients.
	!
	! OUTPUT
	! chebev3d	: value of the approximate function at (x1,x2,x3)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3
	REAL(wp), INTENT(in)	:: x1, x2, x3, coeff( 0:n1, 0:n2, 0:n3), lb1, ub1, lb2, ub2, lb3, ub3
	REAL(wp) :: chebev3d

	! Local variables
	REAL(wp) :: z
	REAL(wp) :: t1( 0:n1 ), t2( 0:n2 ), t3( 0:n3 )
	INTEGER  :: i, j, k


	! 1) Evaluate the Chebyshev polynomials at the points X(1), X(2) and X(3)
	! First dimension
	z = 2.0_wp*( x1 - lb1 ) /( ub1 - lb1 ) - 1.0_wp
	t1( 0 ) = 1.0_wp
	t1( 1 ) = z
	DO i = 2, n1
		t1( i ) = 2.0_wp * z * t1( i - 1 ) - t1( i - 2 )
	END DO

	! Second dimension
	z = 2.0_wp*( x2 - lb2 ) /( ub2 - lb2 ) - 1.0_wp
	t2( 0 ) = 1.0_wp
	t2( 1 ) = z
	DO i = 2, n2
		t2( i ) = 2.0_wp * z * t2( i - 1 ) - t2( i - 2 )
	END DO

	! Third dimension
	z = 2.0_wp*( x3 - lb3 ) /( ub3 - lb3 ) - 1.0_wp
	t3( 0 ) = 1.0_wp
	t3( 1 ) = z
	DO i = 2, n3
		t3( i ) = 2.0_wp * z * t3( i - 1 ) - t3( i - 2 )
	END DO

	! 2) Evaluate the approximation
	chebev3d = 0.0_wp
	DO i = 0, n1
		DO j = 0, n2
			DO k = 0, n3
				chebev3d = chebev3d + coeff( i, j, k ) * t1( i ) * t2( j ) * t3( k )
			END DO
		END DO
	END DO

	END FUNCTION chebev3d


	FUNCTION chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,x3,x4,coeff )

	! PURPOSE
	! To evaluate the 4-dimensional chebyshev approximate function at the values (x1,x2,x3,x4).
	! This function uses Clenshaw algorithm, which is an efficient way of evaluating Chebyshev 
	! polynomials. The one dimensional algorithm is extended to four dimensions (consult my
	! notes).
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! n4	: degree of polynomial dimension 4.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! lb4	: lower bound dimension 4.
	! ub4	: upper bound dimension 4.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! x4	: fourth element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3, 0:n4) with the chebyshev's coefficients.
	!
	! OUTPUT
	! chebev4d	: value of the approximate function at (x1,x2,x3,x4)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3,n4
	REAL(wp), INTENT(in)	:: x1, x2, x3,x4, coeff( 0:n1, 0:n2, 0:n3, 0:n4), lb1, ub1,	lb2,	&
							   & ub2, lb3, ub3, lb4, ub4
	REAL(wp) :: chebev4d

	! Local variables
	REAL(wp) :: z, a( 0:n1, 0:n2, 0:n3 ), b( 0 : n4 + 2 )
	INTEGER  :: i,j,k,l

	! Transform x4 to the interval [-1,1]
	z = 2.0_wp * ( x4 - lb4 )/( ub4 - lb4 ) - 1.0_wp


	! Reduce evaluation to a 2 dimensional problem by summing long third dimension (using 
	! Clenshaw algorithm. )
	DO i = 0, n1
		DO j = 0, n2
			DO k = 0, n3
				b = 0.0_wp
				DO l = n4, 0, -1
					b( l ) = 2.0_wp * z * b( l + 1 ) - b( l + 2 ) + coeff( i, j, k, l )
				END DO
				a( i, j, k ) = b( 0 ) - z * b( 1 )
			END DO
		END DO
	END DO

	! Evaluate chebyshev's approximation as a 3 dimensional problem
	chebev4d = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3, a )

	END FUNCTION chebev4d



	! THE FOLLOWING FUNCTIONS EVALUATE THE DERIVATIVES OF A CHEBYSHEV POLYNOMIAL	


	FUNCTION chebder1d( n, lb, ub, x, coeff )

	! PURPOSE
	! To evaluate the derivative of an n-degree Chebyshev polynomial.
	! Algorithm: the derivative of an n-degree Chebyshev polynomial is an (n-1)-degree
	! Chebyshev polynomial, which coefficients are derived from those of the n-degree
	! polynomial.
	! This subroutine should be used when we need to evaluate the derivative a moderate 
	! number of times. If we need to evaluate the derivative many times, it is convenient 
	! first to CALL the subroutine coeffder1d to find the derived coefficients, and then
	! use the subroutine chebev1d to evaluate the polynomial.
	!
	! INPUTS
	! n		: degree of polynomial
	! lb	: lower bound
	! ub	: upper bound
	! x		: evaluation point of the derivative
	! coeff : array of size (0:n) witht the Chebyshev's coefficients of the original polynomial
	!
	! OUTPUT
	! chebder1d : value of the derivative of the "n" degree chebyshev polynomial

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n
	REAL(wp), INTENT(in)	:: x, coeff( 0 : n ), lb, ub
	REAL(wp) :: chebder1d

	! Local variables
	REAL(wp) :: b( 0 : n + 1 )
	INTEGER ::i

	! Build coefficients of "n-1" degree polynomial from the original "n" degree polynomial
	b(0:n+1) = 0.0_wp			! initialize vector b
	DO i = n-1, 1, -1			! compute recursion
		b( i ) = 2.0_wp * DBLE(i+1) * coeff(i+1) + b(i+2)
	END DO
	b( 0 )  = coeff(1) + b(2)/2.0_wp

	! Normalize to interval [lb,ub]
	b( 0:n-1 ) = ( 2.0_wp / ( ub - lb ) ) * b( 0:n-1 )

	! Evaluate polynomial using Clenshaw recursion
	chebder1d = chebev1d( n - 1 , lb, ub, x, b(0:n-1) )

	END FUNCTION chebder1d


	FUNCTION chebder2d( n1, n2, lb1, ub1, lb2, ub2, x1, x2, coeff, dimder )

	! PURPOSE
	! To evaluate the derivative of a 2 dimensional chebyshev polynomial at a point (x1,x2)
	! Algorithm: let the chebyshev polynomial be of degree n1 in x1, and n2 in x2. Focus 
	! in the partial derivative w.r.t. x. Then this partial derivative can be shown to 
	! be a chebyshev polynomial of degree (n1-1) in x1 and n2 in x2. This is a generalization
	! of the one-dimensional case, and the details can be found on the file chebydoc2.tex.
	! This subroutine should be used when we need to evaluate the derivative a moderate 
	! number of times. If we need to evaluate the derivative many times, it is convenient 
	! first to CALL the subroutine coeffder2d to find the derived coefficients, and then
	! use the subroutine chebev2d to evaluate the polynomial.
	!
	! NOTE: the numerical derivative could be faster.
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = 1, then
	!		  we compute the partial derivative w.r.t. 1st argument. dimder=2 is for the
	!		  second argument.
	!
	! OUTPUT
	! chebder2d	: value of the partial derivative w.r.t. to argument "dimder" at (x1,x2)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, dimder
	REAL(wp), INTENT(in)	:: x1, x2, coeff( 0:n1, 0:n2 ), lb1, ub1, lb2, ub2
	REAL(wp) :: chebder2d

	! Local variables
	REAL(wp) :: b( 0:n1+2, 0:n2+2 )		! working array
	INTEGER ::i, j

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1		
		b = 0.0_wp		! initialize derived coefficients array
		DO j = 0, n2
			DO i = n1-1, 1, -1
				b( i, j ) = 2.0_wp * DBLE(i+1) * coeff( i + 1, j ) + b( i + 2, j )
			END DO
			b( 0, j ) = coeff( 1, j ) + b( 2, j ) / 2.0_wp
		END DO
		! Normalize coefficients to [lb1,ub1]
		b = ( 2.0_wp / ( ub1 - lb1 ) ) * b		
		! Evaluate derived polynomial
		chebder2d = chebev2d( n1-1, n2, lb1, ub1, lb2, ub2, x1, x2, b( 0:n1-1, 0:n2 ) )
	CASE( 2 )						! Partial derivative w.r.t. x2		
		b = 0.0_wp		! initialized derived coefficients array
		DO i = 0, n1
			DO j = n2-1,1, -1
				b( i, j ) = 2.0_wp * DBLE(j+1) * coeff( i, j+1 ) + b( i, j+2 )
			END DO
			b( i, 0 ) = coeff( i, 1 ) + b( i, 2 )/2.0_wp
		END DO
		! Normalize coefficients to [lb2,ub2]
		b = ( 2.0_wp / ( ub2 - lb2 ) ) * b
		! Evaluate derived polynomial
		chebder2d = chebev2d( n1, n2-1, lb1, ub1, lb2, ub2, x1, x2, b( 0:n1, 0:n2-1 ) )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder2d.'
		STOP
	END SELECT

	END FUNCTION chebder2d



	FUNCTION chebder3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3, coeff, dimder )

	! PURPOSE
	! To evaluate the derivative of a 3D chebyshev polynomial at a point (x1,x2,x3)
	! The algorithm is an extension of the one described in the function chebder2d.
	! This subroutine should be used when we need to evaluate the derivative a moderate 
	! number of times. If we need to evaluate the derivative many times, it is convenient 
	! first to CALL the subroutine coeffder3d to find the derived coefficients, and then
	! use the subroutine chebev3d to evaluate the polynomial.
	! Another possibility is to use the numerical derivative function chebder3dn
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = i, then
	!		  we compute the partial derivative w.r.t. ith argument for i=1,2,3
	!
	! OUTPUT
	! chebder3d	: value of the partial derivative w.r.t. to argument "dimder" at (x1,x2,x3)

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, dimder
	REAL(wp), INTENT(in)	:: x1, x2, x3, coeff( 0:n1, 0:n2, 0:n3), lb1, ub1, lb2, ub2, lb3, ub3
	REAL(wp) :: chebder3d

	! Local variables
	REAL(wp) :: b( 0:n1+1, 0:n2+1, 0:n3+1)		! working array
	INTEGER ::i, j, k

	b = 0.0_wp	! initialize working array

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1
		DO k = 0, n3
			DO j = 0, n2
				DO i = n1-1, 1, -1
					b( i, j, k ) = 2.0_wp * DBLE(i+1) * coeff( i + 1, j, k ) + b( i + 2, j, k )
				END DO
				b( 0, j, k ) = coeff( 1, j, k ) + b( 2, j, k ) / 2.0_wp
			END DO
		END DO
			
		! Normalize coefficients to [lb1,ub1]
		b = ( 2.0_wp / ( ub1 - lb1 ) ) * b

		! Evaluate derived polynomial
		chebder3d = chebev3d( n1-1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3,	&
					&		  b( 0:n1-1, 0:n2, 0:n3 ) )
	CASE( 2 )						! Partial derivative w.r.t. x2
		DO k = 0, n3
			DO i = 0, n1
				DO j = n2-1,1, -1
					b( i, j, k ) = 2.0_wp * DBLE(j+1) * coeff( i, j+1, k ) + b( i, j+2, k )
				END DO
				b( i, 0, k ) = coeff( i, 1, k ) + b( i, 2, k )/2.0_wp
			END DO
		END DO

		! Normalize coefficients to [lb2,ub2]
		b = ( 2.0_wp / ( ub2 - lb2 ) ) * b

		! Evaluate derived polynomial
		chebder3d = chebev3d( n1, n2-1, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3,	&
					&		  b( 0:n1, 0:n2-1, 0:n3 ) )
	CASE( 3 )
		DO j = 0, n2
			DO i = 0, n1
				DO k = n3-1,1,-1
					b( i, j, k ) = 2.0_wp * DBLE(k+1) * coeff( i, j, k+1 ) + b( i, j, k+2 )
				END DO
				b( i, j, 0 ) = coeff( i, j, 1 ) + b( i, j, 2 ) / 2.0_wp
			END DO
		END DO

		! Normalize coefficients to [lb3,ub3]
		b = ( 2.0_wp / ( ub3 - lb3 ) ) * b

		! Evaluate derived polynomial
		chebder3d = chebev3d( n1, n2, n3-1, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3,	&
					&		  b( 0:n1, 0:n2, 0:n3-1 ) )

	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder3d.'
		STOP
	END SELECT

	END FUNCTION chebder3d


	FUNCTION chebder4d( n1, n2, n3, n4, lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, x1, x2, x3,x4,	&
					   & coeff, dimder )

	! PURPOSE
	! To evaluate the derivative of a 4-dimensional chebyshev polynomial at a point (x1,x2,x3,x4)
	! Algorithm: let the chebyshev polynomial be of degree n1 in x1, and n2 in x2. Focus 
	! in the partial derivative w.r.t. x. Then this partial derivative can be shown to 
	! be a chebyshev polynomial of degree (n1-1) in x1 and n2 in x2. This is a generalization
	! of the one-dimensional case, and the details can be found on the file chebydoc2.tex.
	! This subroutine should be used when we need to evaluate the derivative a moderate 
	! number of times. If we need to evaluate the derivative many times, it is convenient 
	! first to CALL the subroutine coeffder3d to find the derived coefficients, and then
	! use the subroutine chebev3d to evaluate the polynomial.
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! n4	: degree of polynomial dimension 4.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! lb4	: lower bound dimension 4.
	! ub4	: upper bound dimension 4.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! x4	: fourth element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3, 0:n4) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = i, then
	!		  we compute the partial derivative w.r.t. ith argument for i=1,2,4
	!
	! OUTPUT
	! chebder4d	: value of the partial derivative w.r.t. to argument "dimder" at (x1,x2,x3,x4)


	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, n4, dimder
	REAL(wp), INTENT(in)	:: x1, x2, x3, x4, coeff( 0:n1, 0:n2, 0:n3, 0:n4), lb1, ub1, lb2,	&
							&  ub2,	lb3, ub3, lb4, ub4
	REAL(wp) :: chebder4d

	! Local variables
	REAL(wp) :: b( 0:n1+1, 0:n2+1, 0:n3+1, 0:n4+1)		! working array
	INTEGER ::i, j, k, s

	b = 0.0_wp	! initialize working array

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1
		DO s = 0, n4
			DO k = 0, n3
				DO j = 0, n2
					DO i = n1-1, 1, -1
						b(i,j,k,s) = 2.0_wp * DBLE(i+1) * coeff(i+1,j,k,s) + b(i+2,j,k,s)
					END DO
					b( 0, j, k, s ) = coeff( 1, j, k, s ) + b( 2, j, k, s ) / 2.0_wp
				END DO
			END DO
		END DO
			
		! Normalize coefficients to [lb1,ub1]
		b = ( 2.0_wp / ( ub1 - lb1 ) ) * b

		! Evaluate derived polynomial
		chebder4d = chebev4d( n1-1, n2, n3, n4, lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, x1,		&
							  x2, x3, x4, b( 0:n1-1, 0:n2, 0:n3, 0:n4 ) )
	CASE( 2 )						! Partial derivative w.r.t. x2
		DO s = 0, n4
			DO k = 0, n3
				DO i = 0, n1
					DO j = n2-1,1, -1
						b(i,j,k,s) = 2.0_wp * DBLE(j+1) * coeff(i,j+1,k,s) + b(i,j+2,k,s)
					END DO
					b( i, 0, k, s ) = coeff( i, 1, k, s ) + b( i, 2, k, s )/2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb2,ub2]
		b = ( 2.0_wp / ( ub2 - lb2 ) ) * b

		! Evaluate derived polynomial
		chebder4d = chebev4d( n1, n2-1, n3, n4, lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, x1,		&
							  x2, x3, x4, b( 0:n1, 0:n2-1, 0:n3, 0:n4 ) )
	CASE( 3 )
		DO s = 0, n4
			DO j = 0, n2
				DO i = 0, n1
					DO k = n3-1,1,-1
						b(i,j,k,s) = 2.0_wp * DBLE(k+1) * coeff(i,j,k+1,s) + b(i,j,k+2,s)
					END DO
					b( i, j, 0, s ) = coeff( i, j, 1, s ) + b( i, j, 2, s) / 2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb3,ub3]
		b = ( 2.0_wp / ( ub3 - lb3 ) ) * b

		! Evaluate derived polynomial
		chebder4d = chebev4d( n1, n2, n3-1, n4, lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, x1,		&
							  x2, x3, x4, b( 0:n1, 0:n2, 0:n3-1, 0:n4 ) )
	CASE( 4 )
		DO k = 0, n3
			DO j = 0, n2
				DO i = 0, n1
					DO s = n4-1,1,-1
						b(i,j,k,s) = 2.0_wp * DBLE(s+1) * coeff(i,j,k,s+1) + b(i,j,k,s+2)
					END DO
					b( i, j, k, 0 ) = coeff( i, j, k, 1 ) + b( i, j, k, 2 ) / 2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb4,ub4]
		b = ( 2.0_wp / ( ub4 - lb4 ) ) * b

		! Evaluate derived polynomial
		chebder4d = chebev4d( n1, n2, n3, n4-1, lb1, ub1, lb2, ub2, lb3, ub3, lb4, ub4, x1,		&
							  x2, x3, x4, b( 0:n1, 0:n2, 0:n3, 0:n4-1 ) )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder4d.'
		STOP
	END SELECT

	END FUNCTION chebder4d
	

	SUBROUTINE coeffder1d( n, lb, ub, coeff, dercoeff )
	
	! PURPOSE
	! Given an array coeff of coefficients of an n-degree chebyshev polynomial, this
	! subroutine computes the "( 0: n-1)" coefficients of the derivative of such polynomial.
	!
	! INPUTS
	! n		: degree of polynomial
	! lb	: lower bound
	! ub	: upper bound
	! coeff : array of size (0:n) witht the Chebyshev's coefficients of the original polynomial
	!
	! OUTPUT
	! dercoeff : array of size ( 0:n-1) with the coefficients of the derivative polynomial

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n
	REAL(wp), INTENT(in)	:: coeff( 0 : n ), lb, ub
	REAL(wp), INTENT(out)	:: dercoeff( 0: n-1 )

	! Local variables
	REAL(wp) :: b( 0 : n + 1 )
	INTEGER ::i

	! Build coefficients of "n-1" degree polynomial from the original "n" degree polynomial
	b(0:n+1) = 0.0_wp			! initialize vector b
	DO i = n-1, 1, -1			! compute recursion
		b( i ) = 2.0_wp * DBLE(i+1) * coeff(i+1) + b(i+2)
	END DO
	b( 0 )  = coeff(1) + b(2)/2.0_wp
	
	dercoeff( 0: n-1 ) = ( 2.0_wp / ( ub - lb ) ) * b( 0:n-1 )	

	END SUBROUTINE coeffder1d


	SUBROUTINE coeffder2d( n1, n2, lb1, ub1, lb2, ub2, coeff, dercoeff, dimder )

	! Given an array coeff(0:n1,0:n2) with the coefficients of an (n1,n2)-degree 2 dimensional
	! chebyshev polynoial, this subroutine computes an array dercoeff( 0:n1-1,0:n2 ) or
	! dercoeff( 0:n1, 0:n2-1) with the coefficients of the partial derivative chebyshev
	! polynomial.
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = 1, then
	!		  we compute the partial derivative w.r.t. 1st argument. dimder=2 is for the
	!		  second argument.
	!
	! OUTPUT
	! dercoeff	: array of size either ( 0:n1-1, 0:n2) -if dimder=1- or ( 0:n1, 0:n2-1)
	!			  -if dimder=2- with the coefficients of the derivative polynomial

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, dimder
	REAL(wp), INTENT(in)	:: coeff( 0:n1, 0:n2 ), lb1, ub1, lb2, ub2
	REAL(wp), INTENT(out)  :: dercoeff( :, : )

	! Local variables
	REAL(wp), ALLOCATABLE :: b( :, : )		! working array
	INTEGER ::i, j

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1
		ALLOCATE( b( 0:n1+1, 0:n2 ) )
		b( 0:n1+1, 0:n2 ) = 0.0_wp	! initialize derived coefficients array
		DO j = 0, n2
			DO i = n1-1, 1, -1
				b( i, j ) = 2.0_wp * DBLE(i+1) * coeff( i + 1, j ) + b( i + 2, j )
			END DO
			b( 0, j ) = coeff( 1, j ) + b( 2, j ) / 2.0_wp
		END DO
		! Normalize coefficients to [lb1,ub1]
		b( 0:n1-1, 0:n2 ) = ( 2.0_wp / ( ub1 - lb1 ) ) * b( 0:n1-1, 0:n2 )
		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1-1, 0:n2 )
	CASE( 2 )						! Partial derivative w.r.t. x2
		ALLOCATE( b( 0:n1, 0:n2+1 ) )	
		b( 0:n1, 0:n2+1 ) = 0.0_wp	! initialized derived coefficients array
		DO i = 0, n1
			DO j = n2-1,1, -1
				b( i, j ) = 2.0_wp * DBLE(j+1) * coeff( i, j+1 ) + b( i, j+2 )
			END DO
			b( i, 0 ) = coeff( i, 1 ) + b( i, 2 )/2.0_wp
		END DO
		! Normalize coefficients to [lb2,ub2]
		b( 0:n1, 0:n2-1 ) = ( 2.0_wp / ( ub2 - lb2 ) ) * b( 0:n1, 0:n2-1 )
		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1-1, 0:n2 )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine coeffder2d.'
		STOP
	END SELECT

	END SUBROUTINE coeffder2d


	SUBROUTINE coeffder3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, coeff, dercoeff, dimder )

	! PURPOSE
	! Given an array coeff(0:n1,0:n2, 0:n3) with the coefficients of an (n1,n2,n3)-degree 
	! 3-dimensional chebyshev polynomial, this subroutine computes an array dercoeff of size
	! ( 0:n1-1,0:n2,0:n3 ) or ( 0:n1,0:n2-1,0:n3 ) or ( 0:n1,0:n2,0:n3-1 ) with the 
	! coefficients of the partial derivative chebyshev polynomial.

	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! coeff	: array of size (0:n1, 0:n2, 0:n3) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = i, then
	!		  we compute the partial derivative w.r.t. ith argument for i=1,2,3
	!
	! OUTPUT
	! dercoeff	: array of size (0:n1-1, 0:n2, 0:n3) -if dimder=1-, or ( 0:n1, 0:n2-1, 0:n3)
	!			  -if dimder=2-, or (0:n1, 0:n2, 0:n3-1) with the coefficients of 3D
	!			  derivative polynomial.
	

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, dimder
	REAL(wp), INTENT(in)	:: coeff( 0:n1, 0:n2, 0:n3), lb1, ub1, lb2, ub2, lb3, ub3
	REAL(wp), INTENT(out)  :: dercoeff( :, : , :)

	! Local variables
	REAL(wp) :: b( 0:n1+1, 0:n2+1, 0:n3+1)		! working array
	INTEGER ::i, j, k

	b = 0.0_wp	! initialize working array

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1
		DO k = 0, n3
			DO j = 0, n2
				DO i = n1-1, 1, -1
					b( i, j, k ) = 2.0_wp * DBLE(i+1) * coeff( i + 1, j, k ) + b( i + 2, j, k )
				END DO
				b( 0, j, k ) = coeff( 1, j, k ) + b( 2, j, k ) / 2.0_wp
			END DO
		END DO
			
		! Normalize coefficients to [lb1,ub1]
		b = ( 2.0_wp / ( ub1 - lb1 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1-1, 0:n2, 0:n3 )

	CASE( 2 )						! Partial derivative w.r.t. x2
		DO k = 0, n3
			DO i = 0, n1
				DO j = n2-1,1, -1
					b( i, j, k ) = 2.0_wp * DBLE(j+1) * coeff( i, j+1, k ) + b( i, j+2, k )
				END DO
				b( i, 0, k ) = coeff( i, 1, k ) + b( i, 2, k )/2.0_wp
			END DO
		END DO

		! Normalize coefficients to [lb2,ub2]
		b = ( 2.0_wp / ( ub2 - lb2 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1, 0:n2-1, 0:n3 )

	CASE( 3 )
		DO i = 0, n1
			DO j = 0,n2
				DO k = n3-1,1,-1
					b( i, j, k ) = 2.0_wp * DBLE(k+1) * coeff( i, j, k+1 ) + b( i, j, k+2 )
				END DO
				b( i, j, 0 ) = coeff( i, j, 1 ) + b( i, j, 2 ) / 2.0_wp
			END DO
		END DO

		! Normalize coefficients to [lb3,ub3]
		b = ( 2.0_wp / ( ub3 - lb3 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1, 0:n2, 0:n3-1 )

	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine coeffder3d.'
		STOP
	END SELECT

	END SUBROUTINE coeffder3d



	SUBROUTINE coeffder4d(n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,coeff,dercoeff,dimder )

	! PURPOSE
	! Given an array coeff(0:n1,0:n2, 0:n3, 0:n4) with the coefficients of an (n1,n2,n3,n4)
	! -degree 4D chebyshev polynomial, this subroutine computes an array dercoeff of size
	! ( 0:n1-1,0:n2,0:n3, 0:n4 ) or ( 0:n1,0:n2-1,0:n3, 0:n4 ) or ( 0:n1,0:n2,0:n3-1, 0:n4 ) or
	! ( 0:n1,0:n2,0:n3, 0:n4-1 ) with the coefficients of the partial derivative 
	! chebyshev polynomial.

	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! n4	: degree of polynomial dimension 4.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! lb4	: lower bound dimension 4.
	! ub4	: upper bound dimension 4.
	! coeff	: array of size (0:n1, 0:n2, 0:n3, 0:n4) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = i, then
	!		  we compute the partial derivative w.r.t. ith argument for i=1,2,4
	!
	! OUTPUT
	! dercoeff	: array of size (0:n1-1, 0:n2, 0:n3, 0:n4) -if dimder=1-, or 
	!			 ( 0:n1, 0:n2-1, 0:n3, 0:n4) -if dimder=2-, or (0:n1, 0:n2, 0:n3-1, 0:n4) -if
	!			 dimder=3-, or ( 0:n1, 0:n2, 0:n3, 0:n4-1) -if dimder=4- with the coefficients 
	!			of the 4D derivative polynomial.

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, n4, dimder
	REAL(wp), INTENT(in)	:: coeff(0:n1,0:n2,0:n3,0:n4),lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4
	REAL(wp), INTENT(out)	:: dercoeff(:,:,:,:)


	! Local variables
	REAL(wp) :: b( 0:n1+1, 0:n2+1, 0:n3+1, 0:n4+1)		! working array
	INTEGER ::i, j, k, s

	b = 0.0_wp	! initialize working array

	SELECT CASE( dimder )
	CASE( 1 )						! Partial derivative w.r.t. x1
		DO s = 0, n4
			DO k = 0, n3
				DO j = 0, n2
					DO i = n1-1, 1, -1
						b(i,j,k,s) = 2.0_wp * DBLE(i+1) * coeff(i+1,j,k,s) + b(i+2,j,k,s)
					END DO
					b( 0, j, k, s ) = coeff( 1, j, k, s ) + b( 2, j, k, s ) / 2.0_wp
				END DO
			END DO
		END DO
			
		! Normalize coefficients to [lb1,ub1]
		b = ( 2.0_wp / ( ub1 - lb1 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1-1, 0:n2, 0:n3, 0:n4 )
	CASE( 2 )						! Partial derivative w.r.t. x2
		DO s = 0, n4
			DO k = 0, n3
				DO i = 0, n1
					DO j = n2-1,1, -1
						b(i,j,k,s) = 2.0_wp * DBLE(j+1) * coeff(i,j+1,k,s) + b(i,j+2,k,s)
					END DO
					b( i, 0, k, s ) = coeff( i, 1, k, s ) + b( i, 2, k, s )/2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb2,ub2]
		b = ( 2.0_wp / ( ub2 - lb2 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1, 0:n2-1, 0:n3, 0:n4 )
	CASE( 3 )
		DO s = 0, n4
			DO j = 0, n2
				DO i = 0, n1
					DO k = n3-1,1,-1
						b(i,j,k,s) = 2.0_wp * DBLE(k+1) * coeff(i,j,k+1,s) + b(i,j,k+2,s)
					END DO
					b( i, j, 0, s ) = coeff( i, j, 1, s ) + b( i, j, 2, s) / 2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb3,ub3]
		b = ( 2.0_wp / ( ub3 - lb3 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1, 0:n2, 0:n3-4, 0:n4 )
	CASE( 4 )
		DO k = 0, n3
			DO j = 0, n2
				DO i = 0, n1
					DO s = n4-1,1,-1
						b(i,j,k,s) = 2.0_wp * DBLE(s+1) * coeff(i,j,k,s+1) + b(i,j,k,s+2)
					END DO
					b( i, j, k, 0 ) = coeff( i, j, k, 1 ) + b( i, j, k, 2 ) / 2.0_wp
				END DO
			END DO
		END DO

		! Normalize coefficients to [lb4,ub4]
		b = ( 2.0_wp / ( ub4 - lb4 ) ) * b

		! Construct coefficients of derived polynomial
		dercoeff = b( 0:n1, 0:n2, 0:n3, 0:n4-1 )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine coeffder4d.'
		STOP
	END SELECT

	END SUBROUTINE coeffder4d




	! THE FOLLOWING ROUTINES COMPUTE NUMERICAL DERIVATIVES


	FUNCTION chebder1dn( n, lb, ub, x, coeff )
	
	! PURPOSE
	! To compute the numerical derivative of a 1 dimensional chebyshev polynomial
	! Note: sometimes, the numerical derivative is faster than the routine chebder1d
	!
	! INPUTS
	! n		: degree of chebyshev polynomial
	! lb	: lower bound
	! ub	: upper bound
	! x		: evaluation point of the derivative
	! coeff : array of size (0:n) witht the Chebyshev's coefficients of the original polynomial
	!
	! OUTPUT
	! chebder1dn : value of the numerical derivative of the chebyshev polynomial

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n
	REAL(wp), INTENT(in)	:: x, coeff( 0:n ), lb, ub
	REAL(wp) :: chebder1dn

	!Local variables
	REAL(wp) :: h, hh, xhu, xhl, fu, fl

	h	= max(dabs(x),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
	xhu = x + h
	xhl = x - h
	hh  = xhu - xhl
	fu	= chebev1d( n, lb, ub, xhu, coeff )
	fl	= chebev1d( n, lb, ub, xhl, coeff )

	chebder1dn = ( fu - fl ) / hh

	END FUNCTION chebder1dn


	FUNCTION chebder2dn( n1, n2, lb1, ub1, lb2, ub2, x1, x2, coeff, dimder )

	! PURPOSE
	! To evaluate a numerical partial derivative of a 2 dimensional chebyshev polynomial
	! at the point (x1, x2). 
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. If dimder = 1, then
	!		  we compute the partial derivative w.r.t. 1st argument. dimder=2 is for the
	!		  second argument.
	!
	! OUTPUT
	! chebder2dn	: value of the partial derivative w.r.t. to argument "dimder" at (x1,x2)


	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, dimder
	REAL(wp), INTENT(in)	:: x1, x2, coeff( 0:n1, 0:n2 ), lb1, ub1, lb2, ub2
	REAL(wp) :: chebder2dn

	!Local variables
	REAL(wp) :: h, hh, xhu, xhl, fu, fl
	
	SELECT CASE( dimder )
	CASE(1)
		h = max(dabs(x1),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x1 + h
		xhl = x1 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev2d( n1, n2, lb1, ub1, lb2, ub2, xhu, x2, coeff )
		fl = chebev2d( n1, n2, lb1, ub1, lb2, ub2, xhl, x2, coeff )
	CASE(2)
		h = max(dabs(x2),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x2 + h
		xhl = x2 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev2d( n1, n2, lb1, ub1, lb2, ub2, x1, xhu, coeff )
		fl = chebev2d( n1, n2, lb1, ub1, lb2, ub2, x1, xhl, coeff )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder2dn.'
		STOP
	END SELECT

	chebder2dn = ( fu - fl ) / hh

	END FUNCTION chebder2dn


	FUNCTION chebder3dn( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, x3, coeff, dimder )

	! PURPOSE
	! To evaluate a numerical partial derivative of a 3 dimensional chebyshev polynomial
	! at the point (x1, x2, x3). 
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. ( dimder = 1, first
	!		  argument; dimder = 2, second argument, dimder=3, third argument. )
	!
	! OUTPUT
	! chebder3dn	: value of partial derivative at (x1,x2,x3) with respect to variable dimder

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, dimder
	REAL(wp), INTENT(in)	:: x1, x2, x3, coeff( 0:n1, 0:n2, 0:n3), lb1, ub1, lb2, ub2, lb3, ub3
	REAL(wp) :: chebder3dn

	!Local variables
	REAL(wp) :: h, hh, xhu, xhl, fu, fl
	
	SELECT CASE( dimder )
	CASE(1)
		h = max(dabs(x1),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x1 + h
		xhl = x1 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, xhu, x2, x3, coeff )
		fl = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, xhl, x2, x3, coeff )
	CASE(2)
		h = max(dabs(x2),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x2 + h
		xhl = x2 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, xhu, x3, coeff )
		fl = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, xhl, x3, coeff )
	CASE(3)
		h = max(dabs(x3),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x3 + h
		xhl = x3 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, xhu, coeff )
		fl = chebev3d( n1, n2, n3, lb1, ub1, lb2, ub2, lb3, ub3, x1, x2, xhl, coeff )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder2dn.'
		STOP
	END SELECT

	chebder3dn = ( fu - fl ) / hh

	END FUNCTION chebder3dn



	FUNCTION chebder4dn( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,x3,x4,coeff,dimder )

	! PURPOSE
	! To evaluate a numerical partial derivative of a 3 dimensional chebyshev polynomial
	! at the point (x1, x2, x3). 
	!
	! INPUTS
	! n1	: degree of polynomial dimension 1.
	! n2	: degree of polynomial dimension 2.
	! n3	: degree of polynomial dimension 3.
	! n4	: degree of polynomial dimension 4.
	! lb1	: lower bound dimension 1.
	! ub1	: upper bound dimension 1.
	! lb2	: lower bound dimension 2.
	! ub2	: upper bound dimension 2.
	! lb3	: lower bound dimension 3.
	! ub3	: upper bound dimension 3.
	! lb4	: lower bound dimension 4.
	! ub4	: upper bound dimension 4.
	! x1	: first element of evaluation point.
	! x2	: second element of evaluation point.
	! x3	: third element of evaluation point.
	! x4	: fourth element of evaluation point.
	! coeff	: array of size (0:n1, 0:n2, 0:n3) with the chebyshev's coefficients.
	! dimder: argument the partial derivative is taken with respect to. ( dimder = 1, first
	!		  argument; dimder = 2, second argument, dimder=3, third argument; dimder = 4,
	!		  fourth argument. )
	!
	! OUTPUT
	! chebder4dn	: value of partial derivative at (x1,x2,x3,x4) with respect to 
	!				  variable dimder

	IMPLICIT NONE
	INTEGER, INTENT(in)		:: n1, n2, n3, n4, dimder
	REAL(wp), INTENT(in)	:: x1, x2, x3, x4, coeff( 0:n1, 0:n2, 0:n3), lb1, ub1, lb2, ub2,	&
							   lb3, ub3, lb4, ub4
	REAL(wp) :: chebder4dn

	!Local variables
	REAL(wp) :: h, hh, xhu, xhl, fu, fl
	
	SELECT CASE( dimder )
	CASE(1)
		h = max(dabs(x1),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x1 + h
		xhl = x1 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,xhu,x2,x3,x4,coeff )
		fl = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,xhl,x2,x3,x4,coeff )
	CASE(2)
		h = max(dabs(x2),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x2 + h
		xhl = x2 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,xhu,x3,x4,coeff )
		fl = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,xhl,x3,x4,coeff )
	CASE(3)
		h = max(dabs(x3),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x3 + h
		xhl = x3 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,xhu,x4,coeff )
		fl = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,xhl,x4,coeff )
	CASE(4)
		h = max(dabs(x4),1.0_wp) * EPSILON( 0.0_wp ) ** ( 1.0_wp/3.0_wp )
		xhu = x4 + h
		xhl = x4 - h
		hh	= xhu - xhl		! To make the divided difference representable
		fu = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,x3,xhu,coeff )
		fl = chebev4d( n1,n2,n3,n4,lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4,x1,x2,x3,xhl,coeff )
	CASE DEFAULT
		WRITE(*,*) ' Wrong dimension input in subroutine chebder2dn.'
		STOP
	END SELECT

	chebder4dn = ( fu - fl ) / hh

	END FUNCTION chebder4dn


END MODULE chebyshev_regression