#include "/home/labs/komputerowe_symulacje_numeryczne/numerical_recipes.c/nrutil.h"
#include "/home/labs/komputerowe_symulacje_numeryczne/numerical_recipes.c/nrutil.c"
#include "/home/labs/komputerowe_symulacje_numeryczne/numerical_recipes.c/gaussj.c"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double delta = 0.01;
float xmax=5;
float xmin=-5;
float dx;

double pow2(double x) {
   return x*x;
}

double pow3(double x) {
   return x*x*x;
}

double f1(double x){
   return (1/(1+x*x));
}

double f2(double x){
   return cos(2*x);
}

double fp1(double x){
   return (f1(x+delta) - f1(x-delta)) / (2*delta);
}

double fp2(double x){
   return (f2(x+delta) - f2(x-delta)) / ( 2*delta );
}

double FI(float* x, float X, int i){
   double h = dx;
   double fi = 0;

   if (x[i-2] <= X && X < x[i-1])
      fi = pow3(X - x[i-2]) / pow3(h);

   if (x[i-1] <= X && X < x[i])
      fi = ( pow3(h) + 3.0*pow2(h) * (X-x[i-1]) + 3.0*h * pow2(X-x[i-1]) - 3.0 * pow3(X-x[i-1])) / pow3(h);

   if (x[i] <= X && X < x[i+1])
      fi = ( pow3(h) + 3.0*pow2(h) * (x[i+1]-X) + 3.0*h * pow2(x[i+1]-X) - 3.0 * pow3(x[i+1]-X)) / pow3(h);

   if (x[i+1] <= X && X < x[i+2])
      fi =  pow3(x[i+2] - X) / pow3(h);

   return fi;
}

float s(float X, float* c, int n, float*x){

   double s = 0;
   int i;

   for (i=0; i<=(n+1); i++)
      s += c[i] * FI(x, X, i);
   return s;
}

void fillMatrix(float** m, int n){
   int i,j;
   for (i=1; i<=n; ++i)
      for (j=1; j<=n; ++j)
      {
         m[i][j] = 0;
         if (i==j)
            m[i][j] = 4;
         if (j==(i+1))
            m[i][j] = 1;
         if (j==(i-1))
            m[i][j] = 1;
      }
   m[1][2] = 2;
   m[n][n-1] = 2;
}

void interpolacja(const char* filename, int n, int mode)
{
   dx = 2*xmax/(n-1);
   float** M = matrix(1,n, 1, n);
   fillMatrix(M, n);


   float* xw = vector(-3,n+3);
   float* yw = vector(1,n);

   FILE* fp = fopen(filename, "w");


   int i,j ;
   for(i=-3;i<=(n+3);i++)
      xw[i]=-xmax+dx*(i-1);


   for(i=1;i<=n;i++) {
      if(mode == 0)
        yw[i]=f1(xw[i]);
      else
        yw[i]=f2(xw[i]);
   }


   double alfa, beta;
   if(mode == 0) {
      alfa = fp1(xmin);
      beta = fp1(xmax);
   } else {
      alfa = fp2(xmin);
      beta = fp2(xmax);
   }

   float** W = matrix(1,n,1,1);
   for (i = 1; i<=n; ++i)
      W[i][1] = yw[i];

   W[1][1] = yw[1] + dx*alfa/3;
   W[n][1] = yw[n] - dx*beta/3;

   gaussj(M,n,W,1);
   float* C = vector(0, n+1);
   for (i=1; i<=n; ++i)
      C[i] = W[i][1];

   C[0] = W[2][1] - (dx/3)*alfa;
   C[n+1] = (dx/3)*beta + C[n-1];

   double k;
   for (k=xmin; k<=xmax; k+=delta){
      if(mode == 0)
         fprintf(fp, "%f %f %f\n", k, f1(k), s(k, C, n, xw));
      else
         fprintf(fp, "%f %f %f\n", k, f2(k), s(k, C, n, xw));
   }

}

int main(){
   interpolacja("f1_5.txt",5,0);
   interpolacja("f1_6.txt",6,0);
   interpolacja("f1_10.txt",10,0);
   interpolacja("f1_20.txt",20,0);

   interpolacja("f2_6.txt",6,1);
   interpolacja("f2_7.txt",7,1);
   interpolacja("f2_14.txt",14,1);
}