/* test the fit to a delta function by 
 * Chebyshev and Legendre polynormials */
 
#include <stdio.h>
#include <math.h>
#include "../gkslib/gkslib.h"
#include "../iolib/message.h"
#include "../const/const.h"

#define N 100
 
main()
{
Gwlimit wcsz;
Gwpoint curve[N];
Gwpoint postn;
float *farray1(),**farray2();
float Chebyshev();
float Legpolyn();
float *f,**atb,**ata;
float x0,eta,x;
int i,j,p,pmax;
double add;
char string[120];

  wcsz.xmin=-3.0;
  wcsz.xmax=3.0;
  wcsz.ymin=-2.;
  wcsz.ymax=3.5;
  gstart(&wcsz,1);
  
  gcreateseg(1);
  curve[1].y=curve[0].y=0.0;
  curve[0].x=-1.0;
  curve[1].x=1.0;
  gpolyline(2,curve);
  curve[1].x=curve[0].x=-1.0;
  curve[0].y=-0.5;
  curve[1].y=1.0;
  gpolyline(2,curve);
  gcloseseg();

  printf("Max degree=");
  scanf("%d",&pmax);
  printf("the location of input delta function [0,1]:");
  scanf("%f%*c",&x0);
  printf("damping=");
  scanf("%f%*c",&eta);
  
  gcreateseg(2);
  curve[1].x=curve[0].x=x0;
  curve[0].y=0.0;
  curve[1].y=1.0;
  gpolyline(2,curve);
  
  f=farray1(0,pmax);
  ata=farray2(1,pmax+1,1,pmax+1);      /*RECIPES package*/
  atb=farray2(1,pmax+1,1,1);        /*RECIPES package*/
  
/* Chebyshev */
  
  for(p=0;p<=pmax;p++) f[p]=Chebyshev(p,x0,1);
  
  for(i=1;i<=pmax+1;i++)
  { atb[i][1]=f[i-1];
    for(j=1;j<=pmax+1;j++) ata[i][j]=f[i-1]*f[j-1];
    ata[i][i]+=eta;
  }
  
  gaussj(ata,pmax+1,atb,1);
  
  for(i=0;i<N;i++)
  { curve[i].x=x=-1.0+2.*(F)i/(F)(N-1);
    add=0.0;
    for(p=0;p<=pmax;p++) add+=atb[p+1][1]*Chebyshev(p,x,1);
    curve[i].y=add;
  }
 
  gpolyline(N,curve);
  
/* Legendre */
  
  for(p=0;p<=pmax;p++) f[p]=Legpolyn(p,x0);
  
  for(i=1;i<=pmax+1;i++)
  { atb[i][1]=f[i-1];
    for(j=1;j<=pmax+1;j++) ata[i][j]=f[i-1]*f[j-1];
    ata[i][i]+=eta;
  }
  
  gaussj(ata,pmax+1,atb,1);
  
  for(i=0;i<N;i++)
  { curve[i].x=x=-1.0+2.*(F)i/(F)(N-1);
    add=0.0;
    for(p=0;p<=pmax;p++) add+=atb[p+1][1]*Legpolyn(p,x,1);
    curve[i].y=add;
  }
 
  gsetlinetype(2);
  gpolyline(N,curve);
  
  sprintf(string,"pmax=%d x=%3.1f",pmax,x0);
  postn.x=-0.9;
  postn.y=.9;
  gsetcharheight(0.08);
  gtext(&postn,string);

  free_farray1(f,0);
  free_farray2(ata,1,pmax+1,1);
  free_farray2(atb,1,pmax+1,1);
  
  gcloseseg();
  gwrapup(0);
 
}