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

#include "cosmocalc.h"
#include "geneutils.h"

int main(int narg,char *arg[])
{
 double zmax = 2.5;
 unsigned short flat = 0;
 // -- WMAP
 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
 // -- ouvert matter only
 //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
 // -- plat matter only
 //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.; flat = 1;
 // -- plat lambda only
 //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.; flat = 2;

 // --- la cosmologie
 double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
 cout<<" perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
 CosmoCalc univ(flat,true,zmax);
 univ.SetInteg(perc,dzinc,dzmax,glorder);
 univ.SetDynParam(h100,om0,or0,ol0,w0);
 univ.Print();

 // --- Remplissage de Dloscom en fct de z
 int_4 lmod=10 , nzfin = lmod*100;
 cout<<"Remplissage de loscom vs z: nzfin="<<nzfin<<" , lmod="<<lmod<<endl;
 vector<double> zfin, dlcfin;
 vector<double> zgros, dlcgros;
 for(int_4 i=0;i<=nzfin;i++) {
   double z = i*zmax/(nzfin-1.);
   double loscom = univ.Dloscom(z); // pour test:  loscom = z*(z+4);
   zfin.push_back(z); dlcfin.push_back(loscom);
   if(i%lmod==0) { // ATTENTION il faut etre equidistant en z !
     zgros.push_back(z);
     dlcgros.push_back(loscom);
   }
 }
 nzfin = zfin.size();
 printf("Fin: z(%d) : [%f,%f]  -> loscom : [%f,%f]\n"
       ,nzfin,zfin[0],zfin[nzfin-1],dlcfin[0],dlcfin[nzfin-1]);
 int_4 nzgros = zgros.size();
 printf("Gros: z(%d) : [%f,%f]  -> loscom : [%f,%f]\n"
       ,nzgros,zgros[0],zgros[nzgros-1],dlcgros[0],dlcgros[nzgros-1]);

 // --- Check interpolation en z
 cout<<"Check interpolation en z"<<endl;
 InterpFunc interpgros(zgros[0],zgros[nzgros-1],dlcgros);
 double xnt[10];
 vector<string> noms;
 noms.push_back("z"); noms.push_back("d");
 noms.push_back("d0"); noms.push_back("d1"); noms.push_back("d2");
 NTuple ntg(noms);
 for(int_4 i=0;i<nzfin;i++) {
   unsigned short ok;
   double z = zfin[i];
   xnt[0] = z;
   xnt[1] = dlcfin[i];
   xnt[2] = interpgros(z);
   xnt[3] = interpgros.Linear(z,ok);
   xnt[4] = interpgros.Parab(z,ok);
   ntg.Fill(xnt);
 }

 // --- Check interpolation en los avec inversion
 cout<<"Check interpolation en los avec inversion"<<endl;
 InverseFunc invfun(zgros,dlcgros);
 vector<double> zbackl;
 invfun.ComputeLinear(nzfin,zbackl);
 vector<double> zbackp;
 invfun.ComputeParab(nzfin,zbackp);
 printf("YMin=%f , YMax=%f\n",invfun.YMin(),invfun.YMax());
 printf("Linear ZBack: z(%d) : [%f,%f]\n"
       ,zbackl.size(),zbackl[0],zbackl[zbackl.size()-1]);
 printf("Parab ZBack: z(%d) : [%f,%f]\n"
       ,zbackp.size(),zbackp[0],zbackp[zbackp.size()-1]);

 InterpFunc interpinvl(invfun.YMin(),invfun.YMax(),zbackl);
 InterpFunc interpinvp(invfun.YMin(),invfun.YMax(),zbackp);

 noms.resize(0);
 noms.push_back("z"); noms.push_back("d");
 noms.push_back("zl"); noms.push_back("zp");
 NTuple ntinv(noms);
 for(int_4 i=0;i<nzfin;i++) {
   double z = zfin[i];
   double d = dlcfin[i];
   unsigned short ok;
   xnt[0] = z;
   xnt[1] = d;
   xnt[2] = interpinvl.Parab(d,ok);
   xnt[3] = interpinvp.Parab(d,ok);
   ntinv.Fill(xnt);
 }

 // --- Write PPF
 string tag = "cmvtinterp.ppf";
 cout<<"Write PPF : "<<tag<<endl;
 POutPersist pos(tag);
 tag = "ntg"; pos.PutObject(ntg,tag);
 tag = "ntinv"; pos.PutObject(ntinv,tag);

 return 0;
}

/*
openppf cmvtinterp.ppf

#### Check interpolation
n/plot ntg.d%z ! ! "nsta connectpoints"
n/plot ntg.d0%z ! ! "nsta connectpoints same green"
n/plot ntg.d1%z ! ! "nsta connectpoints same blue"
n/plot ntg.d2%z ! ! "nsta connectpoints same red"

n/plot ntg.d0-d%z ! ! "nsta connectpoints green"
n/plot ntg.d1-d%z ! ! "nsta connectpoints same blue"
n/plot ntg.d2-d%z ! ! "nsta connectpoints same red"
addline 0 0 1000 0

n/plot ntg.d0-d%d ! ! "nsta connectpoints green"
n/plot ntg.d1-d%d ! ! "nsta connectpoints same blue"
n/plot ntg.d2-d%d ! ! "nsta connectpoints same red"
addline 0 0 1000 0

n/plot ntg.(d0-d)/d%d fabs(d)>1e-13 ! "nsta connectpoints green"
n/plot ntg.(d1-d)/d%d fabs(d)>1e-13 ! "nsta connectpoints same blue"
n/plot ntg.(d2-d)/d%d fabs(d)>1e-13 ! "nsta connectpoints same red"
addline 0 0 1000 0

#### Check inversion
n/plot ntinv.d%_nl ! ! "connectpoints"

n/plot ntinv.z%d ! ! "nsta connectpoints"
n/plot ntg.z%d ! ! "nsta connectpoints same red"

n/plot ntinv.z%d ! ! "nsta connectpoints"
n/plot ntinv.zl%d ! ! "nsta connectpoints same blue"
n/plot ntinv.zp%d ! ! "nsta connectpoints same red"

n/plot ntinv.zl-z%d ! ! "nsta connectpoints blue"
n/plot ntinv.zp-z%d ! ! "nsta connectpoints same red"
addline 0 0 10000 0

n/plot ntinv.zl-z%z ! ! "nsta connectpoints blue"
n/plot ntinv.zp-z%z ! ! "nsta connectpoints same red"
addline 0 0 10000 0

n/plot ntinv.(zl-z)/z%z z>0 ! "nsta connectpoints blue"
n/plot ntinv.(zp-z)/z%z z>0 ! "nsta connectpoints same red"
addline 0 0 10000 0

*/