#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 [options] z_redshift"<<endl
     <<" -H h100 -B Ob0 -M Om0 -L Ol0,w0"<<endl
     <<" -k npt,lkmin,lkmax : les valeurs sont en log10"<<endl
     <<" -s scale : on normalise le spectre avec scale"<<endl
     <<" -w : write spectra on ASCII file"<<endl;
}

int main(int narg,char *arg[])
{
 // -- WMAP
 double h100=0.71, Om0=0.267804, Ob0 = 0.0444356, Ol0=0.73,w0=-1.;

 double ns = 1., as = 1.;
 int npt = 10000;
 double lkmin = -3., lkmax=2.;
 double scale = 2.54499e+07;  // normalisation du spectre a z=0 selon SDSS
 bool wrascii = false;

 char c;
  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:")) != -1) {
  switch (c) {
  case 'H' :
    sscanf(optarg,"%lf",&h100);
    break;
  case 'M' :
    sscanf(optarg,"%lf",&Om0);
    break;
  case 'B' :
    sscanf(optarg,"%lf",&Ob0);
    break;
  case 'L' :
    sscanf(optarg,"%lf,%lf",&Ol0,&w0);
    break;
  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 'w' :
    wrascii = true;
    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<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Ol0="<<Ol0<<" w0="<<w0<<endl;
 cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
 cout<<"scale="<<scale<<endl;
 cout<<"zval="<<zval<<endl;

 //--------------------------
 InitialSpectrum pkini(ns,as);
 pkini.SetNorm(scale);

 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.,20.,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 = 100000;
 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);
 }

 //--------------------------
 if(wrascii) {
   cout<<">>>> ASCII"<<endl;
   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),pkz(k,zval),
	     tfnosc2(k),pk0nosc2(k),pkznosc2(k,zval),
	     tfnosc1(k),pk0nosc1(k),pkznosc1(k,zval),
	     tfnob(k),pk0nob(k),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.val%x ! ! "nsta connectpoints"

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

#### Gestion des abscisses
set xk log10(k)
set xk k

#### fct de transfert
zone
n/plot nt.tf%$xk ! ! "nsta connectpoints"
n/plot nt.tfnosc2%$xk ! ! "nsta connectpoints same red"
n/plot nt.tfnosc1%$xk ! ! "nsta connectpoints same blue"
n/plot nt.tfnob%$xk ! ! "nsta connectpoints same green"

n/plot nt.tf/tfnosc2%$xk ! ! "nsta connectpoints red"
n/plot nt.tf/tfnosc1%$xk ! ! "nsta connectpoints same blue"
n/plot nt.tf/tfnob%$xk ! ! "nsta connectpoints same green"
addline -10 1 10 1

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

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

#### Spectre a z
zone
n/plot nt.pk%$xk ! ! "nsta connectpoints"
n/plot nt.pknosc2%$xk ! ! "nsta connectpoints same red"
n/plot nt.pknosc1%$xk ! ! "nsta connectpoints same blue"
n/plot nt.pknob%$xk ! ! "nsta connectpoints same green"

n/plot nt.pk/pknosc2%$xk ! ! "nsta connectpoints red"
n/plot nt.pk/pknosc1%$xk ! ! "nsta connectpoints same blue"
n/plot nt.pk/pknob%$xk ! ! "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%$xk  !  ! "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"

*/