#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <plot.h>

const double PI = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282;
const double F = 2;

double arr[16];

double wave(double t)
{
  double ans = 0, T = 1/F;
  int n;
  for(n = 0;n < 16;++n)
  {
    double tmp;
    if(t == n*T) tmp = 1;
    else tmp = sin(PI*(t-n*T)/T) / (PI*(t-n*T)/T);
    ans += arr[n]*tmp;
  }
  return ans;
}

double original(double t)
{
  //return cos(5 + 1.5)*sin(0.2*PI*t) + 1.3*cos(0.4*PI*t) - 0.9*sin(0.5*PI*t) + 0.5*cos(0.6*PI*t);
  return cos(2*PI*t);
}

void plot(double(*func)(double),double a,double b,double d)
{
  pl_fmove(a,func(a));
  double x;
  for(x = a+d;x <= b;x+=d)
    pl_fcont(x,func(x));
}

int main()
{
  pl_parampl("PAGESIZE","letter");

  int handle = pl_newpl("png",0,stdout,0);
  pl_selectpl(handle);
  pl_openpl();
  pl_fspace(0.0,-2.0,6.0,2.0);
  pl_flinewidth(0.05);
  pl_erase();

  double dt = 0.01;

  pl_pencolorname("red");
  plot(original,0,5,dt);

  int i;
  for(i = 0;i < 16;++i)
    arr[i] = original(i/F);

  pl_pencolorname("blue");
  plot(wave,0,5,dt);

  double t, error = 0;
  for(t = 0;t < 8;t += dt)
  {
    double a = wave(t), b = original(t);
    //printf("a:%f, b:%f, error = %f\n",a,b,fabs(a-b));
    error += fabs(a - b);
  }
  fprintf(stderr,"avg error = %f\n",error/(8/dt));

  pl_closepl();
  pl_selectpl(0);
  pl_deletepl(handle);

  return 0;
}