#include <time.h>
#include <GL/gl.h>
#include <GL/glut.h>
#include <math.h>
#include "splat.h"
#include "filters.c"
#define NX 50
float alphaInt,alphaDeriv;
char normalEstimationMethod[100];

float int_filter[6];
float der_filter[6];
int   int_filter_extent;
int   der_filter_extent;
extern Metrics(int,double *,double *,double *);
char  outfile[200];
/***********************************************************************/
int CreateSmoothDerivative(int EF,int CN,double *data) {
  float filtWeights[6];
  float t;
  int ix,j,lim,extent,scale;
  int nx,numFiltWeights;
  extent=3;
  nx=2*NX*extent+1;
 
  lim=(nx-1)/(2*extent);
  for(ix=0;ix<lim;ix++)     {  
    t=(float)ix/lim;
    numFiltWeights=getFiltWeights((100+10*CN+EF),0,t,filtWeights);
    if(numFiltWeights==6)  {
      for(j=0;j<6;j++)
	data[ix+j*lim]=filtWeights[j];
    }
    else   {
        data[ix]=0;
        for(j=1;j<5;j++)
          data[ix+j*lim]=filtWeights[j-1];
        data[ix+5*lim]=0;
      }
  }        
  t=1.0;
  numFiltWeights=getFiltWeights(100+10*CN+EF,0,t,filtWeights);
    if(numFiltWeights==6)
      data[6*lim]=filtWeights[5];
    else if(numFiltWeights==5)
      data[4*lim]=filtWeights[4];
    return nx;
}  
float smooth_func(int CN, int EF, float x) {
 float filtWeights[6];
 int num;
 float data;
 num=getFiltWeights((100+10*CN+EF),0,x,filtWeights);
 if(num==6) data=filtWeights[3];
 else      data=filtWeights[2];
 return data;
} 
/***********************************************************************/
int CreateCubicKernel(float coeff2,float coeff,double *data)  {  
  int nx;
  float extent;
  
  float a1,b1,c1,a2,b2,c2,scaleInv,x;
  int ix,j,lim,scale;


  extent=2; 
  nx=NX*extent+1; 
  
  scale=(nx-1)/extent;
  scaleInv=1./scale; 
    a1=6+3*coeff+(27.0/6.0)*coeff2;  a2=3*coeff+(1.0/2.0)*coeff2;
    b1=-(2*coeff+6+4*coeff2);        b2=-10*coeff-2*coeff2;
    c1=0;                            c2=8*coeff+2*coeff2;   
    for(ix=0;ix<nx;ix++) {
      x=ix*scaleInv;
      if(x<1.0)
        data[nx+ix-1]=(a1*x+b1)*x;
      else
        data[nx+ix-1]=(a2*x+b2)*x+c2;
       data[nx-ix]=-data[nx+ix-1];
    }
    return nx*2;
}

float cubic_func(float coeff2,float coeff,float x) {
  float a1,b1,c1,a2,b2,c2,scaleInv;
  int ix,j,lim,scale;
  float data;
  a1=6+3*coeff+(27.0/6.0)*coeff2;  a2=3*coeff+(1.0/2.0)*coeff2;
  b1=-(2*coeff+6+4*coeff2);        b2=-10*coeff-2*coeff2;
  c1=0;                            c2=8*coeff+2*coeff2; 
  if(x<1.0)
        data=(a1*x+b1)*x;
      else
        data=(a2*x+b2)*x+c2;
return data;
} 
float Fourier_Cubic(float f) {
  int i,j;
  float I,R,th,data,n;
  R=0;I=0;
  for(n=0;n<2;n+=0.02) {
    th=2*M_PI*f*n;
    R+=(cubic_func(0,-.5,n)*cos(th)-cubic_func(0,-.5,n)*cos(-th));
    I+=(cubic_func(0,-.5,n)*sin(th)-cubic_func(0,-.5,n)*sin(-th));
  }
    data=(double)sqrt(R*R+I*I);
return data;
}
float Fourier_Smooth(float f) {
  int i,j;
  float I,R,th,data,n;
  R=0;I=0;
  for(n=0;n<2;n+=0.02) {
    th=2*M_PI*f*n;
    R+=(smooth_func(0,4,n)*cos(th)-smooth_func(0,4,n)*cos(-th));
    I+=(smooth_func(0,4,n)*sin(th)-smooth_func(0,4,n)*sin(-th));
  }
    data=(double)sqrt(R*R+I*I);
return data;
}/***********************************************************************/
void Draw_Function(void) {
  int i;float scalex;float f,val;
  glColor3f (0.0,1.0,0);
  for(f=0;f<20;f+=.0001) {
    val=Fourier_Smooth(f);
   if( f<.01) printf("\n %f",val);
    glBegin(GL_POINTS);
    glVertex2f((f)*5000/20,val/(2*M_PI));
    glEnd();
  }
  glColor3f (1.0,0.0,0);
  glBegin(GL_LINES);
  glVertex2f(250,-250);     
  glVertex2f(250, 250);
  glEnd();
  glBegin(GL_LINES);
  glVertex2f(0,0);     
  glVertex2f(500, 0);
  glEnd();
  for(i=0;i<=500;i+=125) {
    glBegin(GL_LINES);
    glVertex2f(i,-250);     
    glVertex2f(i, 250);
    glEnd();
  }
 for(i=-250;i<=250;i+=125) {
    glBegin(GL_LINES);
    glVertex2f(0,i);     
    glVertex2f(500, i);
    glEnd();
  }
  glFlush();
}

void Draw_Fourier(double *data, int n) {
  int i;
  for(i=0;i<=n;i++) {
    printf("\n %d %f",i,(float)data[i]);
  }
}

void Draw_Metrics(double *data, int b, int c) {
  int i;
  for(i=0;i<=4;i++) {
    glBegin(GL_POINTS);
    glVertex2f(i,(float)data[i]);
  }
}
void myinit ();
void display ();


void myinit () {

        glClearColor (0, 0, 0, 0);

	glClear (GL_COLOR_BUFFER_BIT);
	glColor3f (1.0,0.0,0);

        glMatrixMode (GL_PROJECTION);
        glLoadIdentity ();
        gluOrtho2D (0, 500 , -25, 25);

        glMatrixMode (GL_MODELVIEW);

}







/***********************************************************************/
main(int argc, char *argv[]) {
  
  FILE *fp1;
  
  glutInit (&argc, argv);
  glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
  glutInitWindowSize (500, 500);

  glutCreateWindow ("GRAPHS");

  glutDisplayFunc (display);
  myinit ();
  glutInitWindowPosition (100,100);
  glutMainLoop ();
  
} 

void display() {
  float C,B;int i;
  float CN,EF;
  
  double *data,*fourier;
  double values[6];int nx;
  data   =(double*)calloc(5000,sizeof(double));
  fourier=(double*)calloc(5000,sizeof(double));

glClear (GL_COLOR_BUFFER_BIT);
 
      nx=NX;
for(i=0;i<=NX/2;i++) data[i]=1;
for(i=0;i<=NX/2;i++) data[NX/2+i]=0;
 

        Draw_Function();

       /*Draw_Fourier(fourier,nx);*/
     /*     Draw_Metrics(values,B,C);*/
 
 
	/*  for (CN=0;CN<=4;CN+=1) {
   for(EF=0;EF<=4;EF+=1) {
     nx=CreateSmoothDerivative(CN,EF,data);
     Fourier(nx,data,fourier);
     Draw_Function(data,nx);*/
	/*     Draw_Metrics(values,B,C);
     }
     }*/

}