#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "srandgen.h"
#include "timing.h"
#include "tvector.h"
#include "ntuple.h"

#include "geneutils.h"

// test vitesse d'interpolation

int main(int narg,char *arg[])
{
  unsigned short typint=0; // 0=near 1=linear 2=parab
  int ntry=1000;
  int npt = 1000;
  double xmin=0., xmax=1.;
  if(narg>1) sscanf(arg[1],"%d",&ntry); 
  if(narg>2) sscanf(arg[2],"%hu",&typint); 
  cout<<"ntry="<<ntry<<" xmin="<<xmin<<" xmax="<<xmax<<" npt="<<npt<<" typint="<<typint<<endl;
  cout<<"total number of calls is: "<<double(npt)*double(ntry)<<endl;

  vector<double> X, Y;
  double dx = (xmax-xmin)/npt;
  double lamb = (xmax-xmin)/20.;
  for(int i=0;i<=npt;i++) {
    double x = xmin + dx*i + dx/3.*drandpm1();
    X.push_back(x);
    double y = sin( 2.*M_PI * x / lamb );
    Y.push_back(y);
  }

  const char *vname[2] = {"x","v"};
  NTuple nt(2,vname);
  double xnt[2];

  InitTim();

  double sumtot = 0.;
  int lptry = ntry/5; if(lptry==0) lptry = 1;
  for(int itry=0;itry<ntry;itry++) {
    bool fill = false;
    if(itry%lptry==0 || itry==ntry-1) fill = true;
    double d = drandpm1()*dx;
    double sum = 0.;
    for(double x=xmin+d;x<xmax+d+dx/5.;x+=dx) {
      double v = InterpTab(x,X,Y,typint);
      sum += v;
      if(fill) {xnt[0] = x; xnt[1] = v; nt.Fill(xnt);}
    }
    sumtot += sum/(double)npt;
  }
  cout<<"sum="<<sumtot/(double)ntry<<endl;

  PrtTim("Fin du job");

  {
  POutPersist pos("cmvtinttab.ppf");
  TVector<r_8> VX(X);
  TVector<r_8> VY(Y);
  pos.PutObject(VX,"x");
  pos.PutObject(VY,"y");
  pos.PutObject(nt,"nt");
  }


 return 0;
}

/*
openppf cmvtinttab.ppf

n/plot x.val%n ! ! "plusmarker3"
n/plot y.val%n ! ! "plusmarker3"
vecplot x y "plusmarker3"

n/plot nt.v%x ! ! "plusmarker3"
vecplot x y "circlemarker5 same red"

 */