/*******************************************************
*    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].                              *
*                                                      *
*                C++ Release 1.1 By J-P Moreau, Paris. *
*                        (www.jpmoreau.fr)             *
* ---------------------------------------------------- *
* Release 1.1: added functions CHINT and CHDER.        *
*******************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  51
#define  HALF  0.5
#define  ONE   1.0
#define  QUART 0.25
#define  TWO   2.0
#define  ZERO  0.0
#define  PI    4.0*atan(1.0)

double FUNC(double x);

void CHEBFT(double A, double B, double *C, int N)  {
/*******************************************************
* 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.                  *
*******************************************************/
  double SUM, F[NMAX];
  double BMA,BPA,FAC, Y;
  int J,K;
  BMA=HALF*(B-A); BPA=HALF*(B+A);
  for (K=1; K<=N; K++) {
    Y=cos(PI*(K-HALF)/N);
    F[K]=FUNC(Y*BMA+BPA);
  }
  FAC=TWO/N;
  for (J=1; J<=N; J++) {
    SUM=ZERO;
    for (K=1; K<=N; K++)
      SUM += F[K]*cos((PI*(J-1))*((K-HALF)/N));
    C[J]=FAC*SUM;
  }
}

double CHEBEV(double A, double B, double *C, int M, double 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.                                    *
*********************************************************/
  double D,DD,SV,Y,Y2;
  int J;
  if ((X-A)*(X-B) > ZERO) printf("\n X not in range.\n\n");
  D=ZERO; DD=ZERO;
  Y=(TWO*X-A-B)/(B-A);  // change of variable
  Y2=TWO*Y;
  for (J=M; J>1; J--) {
    SV=D;
    D=Y2*D-DD+C[J];
    DD=SV;
  }
  return (Y*D-DD+HALF*C[1]);
}

void CHINT(double A,double B, double *C, double *CINT, int 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.             *
*********************************************************/
  double CON,FAC,SUM;
  int J;
  CON=QUART*(B-A);
  SUM=ZERO;
  FAC=ONE;
  for (J=2; J<N; J++) {
    CINT[J]=CON*(C[J-1]-C[J+1])/(J-1);
    SUM += FAC*CINT[J];
    FAC=-FAC;
  }
  CINT[N]=CON*C[N-1]/(N-1);
  SUM += FAC*CINT[N];
  CINT[1]=TWO*SUM;      //set the constant of integration
}

void CHDER(double A,double B, double *C, double *CDER, int 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.                              *
*********************************************************/
  double CON;
  int J;
  CDER[N]=ZERO;
  CDER[N-1]=TWO*(N-1)*C[N];
  if (N >= 3)
    for (J=N-2; J>0; J--)
      CDER[J]=CDER[J+2]+TWO*J*C[J+1];
  CON=TWO/(B-A);
  for (J=1; J<=N; J++)      //normalize to interval B - A
    CDER[J] *= CON;
}

// end of file Chebyshe.cpp