#include <iostream>
#include <cmath>
#include "smesh.h"
#include "sfunction.h"

using namespace std;

double func1(double x)
{  return exp(-fabs(x))/sqrt(fabs(x));}

double func2(double x)
{  return sin(x)*exp(-x);}

double func3(double x)
{  return sin(1/x);}

void CheckSplineError(const mesh1D& om, const function1D<double>& f,
		      const mesh1D& x, const function1D<double>& g, double fp0, double fpN)
{
  spline1D<double> spln(om, f, fp0, fpN);
  clog.precision(16);
  clog<<"# Integral of splined function is "<<spln.integrate()<<endl;

  tint position = om.InitInterpLeft();
  for (int i=0; i<x.size(); i++){
    cout<<x[i]<<" "<<g[i]-spln(om.InterpLeft(x[i],position))<<" "<<g[i]<<endl;
  }
}

void FourierTranformSpline(const mesh1D& spline_x, const function1D<double>& spline_f,
			   const mesh1D& x, double fp0, double fpN)
{
  spline1D<double> spln(spline_x, spline_f, fp0, fpN);
  clog.precision(16);
  clog<<"# Integral of splined function is "<<spln.integrate()<<endl;
  
  tint position = spline_x.InitInterpLeft();
  for (int i=0; i<x.size(); i++){
    dcomplex four = spln.Fourier(x[i], spline_x);
    cout<<x[i]<<" "<<four.real()<<" "<<four.imag()<<endl;
  }
}

int main(int argc, char *argv[], char *env[])
{
  double x0=1e-6, x1=1, x2=10; // mesh parameters
  double fp0 = -5e8, fpN = -1.5e-5; // derivative at zero and N
  int N0 = 200;
  int N1 = 2000;
  int task=0;
  
  int i=0;
  while (++i<argc){
    std::string str(argv[i]);
    if (str=="-task" && i<argc-1) task = atoi(argv[++i]);
    if (str=="-x0" && i<argc-1)   x0   = atof(argv[++i]);
    if (str=="-x1" && i<argc-1)   x1   = atof(argv[++i]);
    if (str=="-x2" && i<argc-1)   x2   = atof(argv[++i]);
    if (str=="-N0" && i<argc-1)   N0   = atoi(argv[++i]);
    if (str=="-N1" && i<argc-1)   N1   = atoi(argv[++i]);
    if (str=="-fp0" && i<argc-1)  fp0  = atof(argv[++i]);
    if (str=="-fpN" && i<argc-1)  fpN  = atof(argv[++i]);
    if (str=="-h" || str=="--help"){
      std::clog<<"**************** Spline test program *********\n";
      std::clog<<"**   Copyright Kristjan Haule, 24.1.2006    **\n";
      std::clog<<"**********************************************\n";
      std::clog<<"Usage: spline [Options]\n";
      std::clog<<"Options:   -task      Check error: 0, Compute Fourier: 1 ("<<task<<")\n";
      std::clog<<"           -x0        start of log mesh  ("<<x0<<")\n";
      std::clog<<"           -x1        start of tan mesh  ("<<x1<<")\n";
      std::clog<<"           -x2        end of mesh  ("<<x2<<")\n";
      std::clog<<"           -N0        Number of points in the mesh  ("<<N0<<")\n";
      std::clog<<"           -N1        Number of points in the mesh used for checking error of the interpolation ("<<N1<<")\n";
      std::clog<<"           -fp0       Derivative at first point ("<<fp0<<")\n";
      std::clog<<"           -fpN       Derivative at last point ("<<fpN<<")\n";
      return 0;
    }
  }
  
  mesh1D om;
  om.MakePositiveLogTanMesh(N0, x0, x1, x2); // Grid of 200 points is constructed
  
  function1D<double> f(om.size());
  for (int i=0; i<om.size(); i++) f[i] = func3(om[i]); // Function evaluated on this mesh
  

  mesh1D x;
  x.MakePositiveLogTanMesh(N1, x0, x1, x2); // 10-times bigger grid is created
  function1D<double> g(x.size());
  for (int i=0; i<x.size(); i++) g[i] = func3(x[i]); // Function also evaluated on the larger grid

  switch (task){
  case 0 : {CheckSplineError(om, f, x, g, fp0, fpN); break;}
  case 1 : {FourierTranformSpline(om, f, om, fp0, fpN); break;}
  }
  
  return 0;
  
}