#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 "histos.h"
#include "tvector.h"
#include "ntuple.h"
#include "perandom.h"

#include "constcosmo.h"
#include "pkspectrum.h"
#include "geneutils.h"

void usage(void);
void usage(void) {
 cout<<"cmvtstpk -k npt,lkmin,lkmax -s scale zval"<<endl
     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
     <<" -s scale : on multiplie pkz par scale"<<endl;
}

int main(int narg,char *arg[])
{
 double ob0 = 0.0444356;
 // -- 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.;

 double ns = 1., as = 1.;

 int npt = 1000;
 double lkmin = -4., lkmax=2.;
 double scale = 1.;

 char c;
  while((c = getopt(narg,arg,"hk:s:")) != -1) {
  switch (c) {
  case 'k' :
    sscanf(optarg,"%d,%lf,%lf",&npt,&lkmin,&lkmax);
    if(npt<=0) npt=1000;
    if(lkmax<lkmin) {lkmin=-4.; lkmax=2.;}
    break;
  case 's' :
    sscanf(optarg,"%lf",&scale);
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }
 double dlk=(lkmax-lkmin)/npt;

 // Fill z value into list
 double zval = 0.;
 if(optind<narg) zval = atof(arg[optind]);

 cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
 cout<<"scale="<<scale<<endl;
 cout<<"zval="<<zval<<endl;

 //--------------------------
 InitialSpectrum pkini(ns,as);

 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
 cout<<"kpeak="<<tf.KPeak()<<endl;
 TransfertEisenstein tfnosc2(tf); tfnosc2.SetNoOscEnv(2);
 TransfertEisenstein tfnosc1(tf); tfnosc1.SetNoOscEnv(1);
 TransfertEisenstein tfnob(h100,om0,0.,T_CMB_Par,true);

 GrowthFactor d1(om0,ol0);
 cout<<"GrowthFactor: "<<d1(zval)<<endl;

 PkSpectrum0 pk0(pkini,tf);
 PkSpectrum0 pk0nosc2(pkini,tfnosc2);
 PkSpectrum0 pk0nosc1(pkini,tfnosc1);
 PkSpectrum0 pk0nob(pkini,tfnob);

 PkSpectrumZ pkz(pk0,d1,zval);
 PkSpectrumZ pkznosc2(pk0nosc2,d1,zval);
 PkSpectrumZ pkznosc1(pk0nosc1,d1,zval);
 PkSpectrumZ pkznob(pk0nob,d1,zval);

 //--------------------------
 Histo hd1(0.,1000.,10000); hd1.ReCenterBin();
 FuncToHisto(d1,hd1,false);

 Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
 FuncToHisto(pkz,hpkz,true);
 TVector<r_8> vpkz(npt);
 FuncToVec(pkz,vpkz,lkmin,lkmax,true);

 FunRan talea(hpkz,true);
 Histo halea(hpkz); halea.Zero();
 int nalea = 1000000;
 for(int i=0;i<nalea;i++) halea.Add(talea.Random());
 halea *= hpkz.Sum()/halea.Sum();

 //--------------------------
 const int n = 14;
 char *vname[n] = {
   "k","pkini",
   "tf","pk0","pk",
   "tfnosc2","pk0nosc2","pknosc2",
   "tfnosc1","pk0nosc1","pknosc1",
   "tfnob","pk0nob","pknob"
 };
 NTuple nt(n,vname);
 double xnt[n];

 for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
   double k = pow(10.,lk);
   xnt[0] = k;
   xnt[1] = pkini(k);
   xnt[2] = tf(k);
   xnt[3] = pk0(k);
   xnt[4] = pkz(k,zval);
   xnt[5] = tfnosc2(k);
   xnt[6] = pk0nosc2(k);
   xnt[7] = pkznosc2(k,zval);
   xnt[8] = tfnosc1(k);
   xnt[9] = pk0nosc1(k);
   xnt[10] = pkznosc1(k,zval);
   xnt[11] = tfnob(k);
   xnt[12] = pk0nob(k);
   xnt[13] = pkznob(k,zval);
   nt.Fill(xnt);
 }

 //--------------------------
 cout<<">>>> ASCII"<<endl;
 if(1) {
   FILE * fdata = fopen("cmvtstpk.data","w");
   fprintf(fdata,"# z= %g : k(Mpc^-1) pkini, tf pk0 pkz, tfnosc2 pk0nosc2 pkznosc2, ",zval);
   fprintf(fdata,"tfnosc1 pk0nosc1 pkznosc1, tfnob pk0nob pkznob\n");
   for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
     double k = pow(10.,lk);
     fprintf(fdata,"%e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
	     k,pkini(k),
	     tf(k),pk0(k),scale*pkz(k,zval),
	     tfnosc2(k),pk0nosc2(k),scale*pkznosc2(k,zval),
	     tfnosc1(k),pk0nosc1(k),scale*pkznosc1(k,zval),
	     tfnob(k),pk0nob(k),scale*pkznob(k,zval));
   }
   fclose(fdata);
 }

 //--------------------------
 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvtstpk.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
 tag = "hd1"; pos.PutObject(hd1,tag);
 tag = "hpkz"; pos.PutObject(hpkz,tag);
 tag = "vpkz"; pos.PutObject(vpkz,tag);
 tag = "halea"; pos.PutObject(halea,tag);
 return 0;
}

/*
openppf cmvtstpk.ppf

#### growth-factor
zone
n/plot hd1.log10(val)%log10(x) val>0.&&x>0. ! "nsta connectpoints"

#### Spectre initial
zone
n/plot nt.log10(pkini)%log10(k) pkini>0. ! "nsta crossmarker3"

#### fct de transfert
zone
n/plot nt.tf%log10(k) ! ! "nsta connectpoints"
n/plot nt.tfnosc2%log10(k) ! ! "nsta connectpoints same red"
n/plot nt.tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
n/plot nt.tfnob%log10(k) ! ! "nsta connectpoints same green"

n/plot nt.tf/tfnosc2%log10(k) ! ! "nsta connectpoints red"
n/plot nt.tf/tfnosc1%log10(k) ! ! "nsta connectpoints same blue"
n/plot nt.tf/tfnob%log10(k) ! ! "nsta connectpoints same red"
addline -10 1 10 1

#### Spectre a z=0
zone
n/plot nt.pk0%log10(k) ! ! "nsta connectpoints"
n/plot nt.pk0nosc2%log10(k) ! ! "nsta connectpoints same red"
n/plot nt.pk0nosc1%log10(k) ! ! "nsta connectpoints same blue"
n/plot nt.pk0nob%log10(k) ! ! "nsta connectpoints same green"

# Check
zone 2 2
n/plot nt.pk0/pkini-tf*tf%log10(k) pkini>0 ! "nsta crossmarker3"
n/plot nt.pk0nosc2/pkini-tfnosc2*tfnosc2%log10(k) pkini>0 ! "nsta crossmarker3"
n/plot nt.pk0nosc1/pkini-tfnosc1*tfnosc1%log10(k) pkini>0 ! "nsta crossmarker3"
n/plot nt.pk0nob/pkini-tfnob*tfnob%log10(k) pkini>0 ! "nsta crossmarker3"

#### Spectre a z
zone
n/plot nt.pk%log10(k) ! ! "nsta connectpoints"
n/plot nt.pknosc2%log10(k) ! ! "nsta connectpoints same red"
n/plot nt.pknosc1%log10(k) ! ! "nsta connectpoints same blue"
n/plot nt.pknob%log10(k) ! ! "nsta connectpoints same green"

n/plot nt.pk/pknosc2%log10(k) ! ! "nsta connectpoints red"
n/plot nt.pk/pknosc1%log10(k) ! ! "nsta connectpoints same blue"
n/plot nt.pk/pknob%log10(k) ! ! "nsta connectpoints same green"
addline -10 1 10 1

#### Le spectre version Delta^2
set D2 k*k*k*pk/(2*M_PI*M_PI)
n/plot nt.$D2%log10(k)  !  ! "nsta crossmarker3 connectpoints"

#### Test des transferts dans Histo et TVector
zone 1 2
n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"
disp hpkz "same red"

disp vpkz
c++exec cout<<hpkz(0)<<" =?= "<<vpkz(0)<<endl;

zone
disp hpkz
disp halea "same red"

*/