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

using namespace std;

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

int main()
{
  double x0=1e-6, x1=1, x2=10; // mesh parameters
  
  int N0 = 200;
  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] = func(om[i]); // Function evaluated on this mesh

  int N1 = 2000;
  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] = func(x[i]); // Function also evaluated on the large grid
  

  tint position = om.InitInterpLeft();   // Initialization of interpolation with increasing argument
  for (int i=0; i<x.size(); i++){
    intpar p = om.InterpLeft(x[i],position); // calculates location in ordered table and interpolating coefficient
    cout<<x[i]<<" "<<f(p)-g[i]<<" "<<g[i]<<endl; // operator()(intpar p) actually interpolates function f
  }
  return 0;
}