#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
<<"  -F filename : read also CAMB 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;
 string fcmbfile = "";

 char c;
  while((c = getopt(narg,arg,"hwk:s:H:M:B:L:F:")) != -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 'F' :
    fcmbfile = optarg;
    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<<" Ocdm="<<Om0-Ob0<<" Ol0="<<Ol0<<" w0="<<w0<<endl;
 cout<<"lkmin="<<lkmin<<" lkmax="<<lkmax<<" npt="<<npt<<" dlk="<<dlk<<endl;
 cout<<"scale="<<scale<<endl;
 cout<<"zval="<<zval<<endl;
 cout<<"fcmbfile="<<fcmbfile<<endl;

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

 TransfertEisenstein tf(h100,Om0-Ob0,Ob0,T_CMB_Par,false);
   double tf0 = tf(0);
   cout<<"kpeak="<<tf.KPeak()<<endl;
 TransfertEisenstein tfnosc2(tf); tfnosc2.SetNoOscEnv(2);
   double tfnosc20 = tfnosc2(0);
 TransfertEisenstein tfnosc1(tf); tfnosc1.SetNoOscEnv(1);
   double tfnosc10 = tfnosc1(0);
 TransfertEisenstein tfnob(h100,Om0,0.,T_CMB_Par,true);
   double tfnob0 = tfnob(0);
 cout<<"tf(0)="<<tf0<<" tfnosc2(0)="<<tfnosc20<<" tfnosc1(0)="<<tfnosc10<<" tfnob(0)="<<tfnob0<<endl;

 GrowthEisenstein d1(Om0,Ol0);
 cout<<"GrowthFactor: D1(0)="<<d1(0.)<<" D1("<<zval<<") = "<<d1(zval)<<endl;

 PkSpecCalc pk0(pkini,tf,d1,0.);
 PkSpecCalc pk0nosc2(pkini,tfnosc2,d1,0.);
 PkSpecCalc pk0nosc1(pkini,tfnosc1,d1,0.);
 PkSpecCalc pk0nob(pkini,tfnob,d1,0.);

 PkSpecCalc pkz(pkini,tf,d1,zval);
 PkSpecCalc pkznosc2(pkini,tfnosc2,d1,zval);
 PkSpecCalc pkznosc1(pkini,tfnosc1,d1,zval);
 PkSpecCalc pkznob(pkini,tfnob,d1,zval);

 PkTabulate* pkcamb = NULL;
 if(fcmbfile.size()>0) {
   cout<<endl<<"PkTabulate for CAMB"<<endl;
   pkcamb = new PkTabulate;
   pkcamb->ReadCAMB(fcmbfile,h100,0.);
   pkcamb->SetInterpTyp(2);
   if(pkcamb->NPoints()==0) {delete pkcamb; pkcamb = NULL;}
 }

 //--------------------------
 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 = 15;
 const char *vname[n] = {
   "k","pkini",
   "tf","pk0","pk",
   "tfnosc2","pk0nosc2","pknosc2",
   "tfnosc1","pk0nosc1","pknosc1",
   "tfnob","pk0nob","pknob","pktab"
 };
 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);
   if(pkcamb != NULL) xnt[14] = (*pkcamb)(k); else xnt[14] = 0;
   nt.Fill(xnt);
 }

   if(pkcamb != NULL) {delete pkcamb; pkcamb = NULL;}

 //--------------------------
 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");
   int n = 0;
   for(double lk=lkmin;lk<lkmax+dlk/2.;lk+=dlk) {
     double k = pow(10.,lk);
     fprintf(fdata,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
	     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));
     n++;
   }
   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"
n/plot hd1.1./(1.+x)%x ! ! "nsta connectpoints same red"

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

#### fct de transfert
zone
n/plot nt.tf%k k>0&&tf>0 ! "nsta connectpoints logx logy"
n/plot nt.tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints same red logx logy"
n/plot nt.tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx logy"
n/plot nt.tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx logy"

n/plot nt.tf/tfnosc2%k k>0&&tfnosc2>0 ! "nsta connectpoints red logx"
n/plot nt.tf/tfnosc1%k k>0&&tfnosc1>0 ! "nsta connectpoints same blue logx"
n/plot nt.tf/tfnob%k k>0&&tfnob>0 ! "nsta connectpoints same green logx"
addline -10 1 10 1

#### Spectre a z=0
zone
n/plot nt.pk0%k k>0&&pk0>0 ! "nsta connectpoints logx logy"
n/plot nt.pk0nosc2%k k>0&&pk0nosc2>0 ! "nsta connectpoints same red logx logy"
n/plot nt.pk0nosc1%k k>0&&pk0nosc1>0 ! "nsta connectpoints same blue logx logy"
n/plot nt.pk0nob%k k>0&&pk0nob>0 ! "nsta connectpoints same green logx logy"

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

#### Spectre a z
zone
n/plot nt.pk%k k>0&&pk>0 ! "nsta connectpoints logx logy"
n/plot nt.pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx logy"
n/plot nt.pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx logy"
n/plot nt.pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx logy"

n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints red logx"
n/plot nt.pk/pknosc1%k k>0&&pknosc1>0 ! "nsta connectpoints same blue logx"
n/plot nt.pk/pknob%k k>0&&pknob>0 ! "nsta connectpoints same green logx"
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%k  k>0  ! "nsta crossmarker3 connectpoints logx"

#### 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"

#### Le spectre tabule
n/plot nt.pk%k pk>0&&k>0 ! "nsta connectpoints logx logy"
n/plot nt.pktab%k pktab>0&&k>0 ! "nsta connectpoints same red logx logy"

n/plot nt.(pktab-pk)%k k>0&&pk>0 ! "nsta connectpoints"
n/plot nt.(pktab-pk)/pk%k k>0&&pk>0 ! "nsta connectpoints"

n/plot nt.pk/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints logx"
n/plot nt.pktab/pknosc2%k k>0&&pknosc2>0 ! "nsta connectpoints same red logx"

*/