// Check cmv / Wayne Hu coding of the transfert functions
// ATTENTION: ne pas appeler hu_power.c ici car il ecrase les initialisations de hu_tf_fit.c
#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "ntuple.h"

#include "constcosmo.h"
#include "pkspectrum.h"
#include "hu_tf_fit.h"

int main(int narg,char *arg[])
{
 double k1 = 1e-4, k2 = 10.; int_4 npt = 10000;
 double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
 //double h100 = 0.71, Om0 = 0.2736, Ob0 = 0.0456, tcmb = T_CMB_Par; //WMAP5+BAO+SN
 //double h100 = 0.5, Om0 = 0.2, Ob0 = 0.1, tcmb = T_CMB_Par; //papier de Hu
 //double h100 = 0.71, Om0 = 1.0, Ob0 = 0.0444356, tcmb = T_CMB_Par; // Om=1 Ob<<1

 cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
 cout<<"Om0*h^2="<<Om0*h100*h100<<" Ob0*h^2="<<Ob0*h100*h100<<endl;
 double fbaryon = Ob0/Om0;
 cout<<"fbaryon="<<fbaryon<<endl;

 cout<<endl<<"TFset_parameters"<<endl;
 TFset_parameters(Om0*h100*h100,fbaryon,tcmb);
 TFprint_parameters();
 double soundfithu = TFsound_horizon_fit(Om0,fbaryon,h100)/h100;
 cout<<"TFsound_horizon_fit = "<<soundfithu<<" Mpc"<<endl;

 cout<<endl<<"TransfertEisenstein with baryon"<<endl;
 TransfertEisenstein Tfull(h100,Om0-Ob0,Ob0,tcmb,false,1);

 cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.2"<<endl;
 TransfertEisenstein Tnowig(h100,Om0-Ob0,Ob0,tcmb,false,1);
 Tnowig.SetNoOscEnv(2);

 cout<<endl<<"TransfertEisenstein without baryon"<<endl;
 TransfertEisenstein Tnob(h100,Om0,0.,tcmb,true,1);

 cout<<endl<<"Compare kpeak"<<endl;
 double kpeak = Tfull.KPeak();
 double kpeakhu = TFk_peak(Om0,fbaryon,h100)*h100;
 cout<<"kpeak = "<<kpeak<<" "<<kpeakhu<<"  Mpc^-1  -> diff "<<kpeak-kpeakhu<<endl;

 const int n = 11;
 const char *vname[n] = {"k","tf","tfnw","tfnb","tfc","tfb"
                  ,"hutf","hutfnw","hutfnb","hutfc","hutfb"};
 NTuple nt(n,vname);
 double xnt[n];

 double kmax[5]={0,0,0,0,0};
 double tmax[5]={0.,0.,0.,0.,0.};
 double lnk1 = log10(k1), lnk2=log10(k2), dlnk=(lnk2-lnk1)/npt;
 for(double lnk=lnk1;lnk<lnk2+dlnk/2.;lnk+=dlnk) {
   double k = pow(10.,lnk);
   //if(fabs(lnk-log10(0.0186209))>1e-5) continue;
   Tfull.SetReturnPart(TransfertEisenstein::ALL); double tf = Tfull(k);
   double tfnw = Tnowig(k);
   double tfnb = Tnob(k);
   Tfull.SetReturnPart(TransfertEisenstein::CDM); double tfc = Tfull(k);
   Tfull.SetReturnPart(TransfertEisenstein::BARYON); double tfb = Tfull(k);
   double hutfc, hutfb;
   double hutf = TFfit_onek(k,&hutfb,&hutfc); 
   double hutfnw = TFnowiggles(Om0,fbaryon,h100,tcmb,k/h100);
   double hutfnb = TFzerobaryon(Om0,h100,tcmb,k/h100);


   if(fabs(tf-hutf)>tmax[0])     {kmax[0]=k; tmax[0]=fabs(tf-hutf);}
   if(fabs(tfnw-hutfnw)>tmax[1]) {kmax[1]=k; tmax[1]=fabs(tfnw-hutfnw);}
   if(fabs(tfnb-hutfnb)>tmax[2]) {kmax[2]=k; tmax[2]=fabs(tfnb-hutfnb);}
   if(fabs(tfc-hutfc)>tmax[3])   {kmax[3]=k; tmax[3]=fabs(tfc-hutfc);}
   if(fabs(tfb-hutfb)>tmax[4])   {kmax[4]=k; tmax[4]=fabs(tfb-hutfb);}

   xnt[0]=k;
   xnt[1]=tf; xnt[2]=tfnw; xnt[3]=tfnb; xnt[4]=tfc; xnt[5]=tfb;
   xnt[6]=hutf; xnt[7]=hutfnw; xnt[8]=hutfnb; xnt[9]=hutfc; xnt[10]=hutfb;
   nt.Fill(xnt);
 }

 cout<<"Maximum difference for:"<<endl;
 cout<<"tf:   "<<tmax[0]<<"  for k="<<kmax[0]<<endl;
 cout<<"tfnw: "<<tmax[1]<<"  for k="<<kmax[1]<<endl;
 cout<<"tfnb: "<<tmax[2]<<"  for k="<<kmax[2]<<endl;
 cout<<"tfc:  "<<tmax[3]<<"  for k="<<kmax[3]<<endl;
 cout<<"tfb:  "<<tmax[4]<<"  for k="<<kmax[4]<<endl;

 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvchkwhu.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);

 return 0;
}

/*
openppf cmvchkwhu.ppf

set k k
set k log10(k)

######################
zone 3 1
n/plot nt.tf%$k ! ! "nsta connectpoints"
n/plot nt.hutf%$k ! ! "nsta connectpoints red same"
n/plot nt.tfc%$k ! ! "nsta connectpoints"
n/plot nt.hutfc%$k ! ! "nsta connectpoints red same"
n/plot nt.tfb%$k ! ! "nsta connectpoints"
n/plot nt.hutfb%$k ! ! "nsta connectpoints red same"
zone 2 1
n/plot nt.tfnw%$k ! ! "nsta connectpoints"
n/plot nt.hutfnw%$k ! ! "nsta connectpoints red same"
n/plot nt.tfnb%$k ! ! "nsta connectpoints"
n/plot nt.hutfnb%$k ! ! "nsta connectpoints red same"

######################
zone 3 1
n/plot nt.tf-hutf%$k ! ! "nsta connectpoints"
n/plot nt.tfc-hutfc%$k ! ! "nsta connectpoints"
n/plot nt.tfb-hutfb%$k ! ! "nsta connectpoints"
zone 2 1
n/plot nt.tfnw-hutfnw%$k ! ! "nsta connectpoints"
n/plot nt.tfnb-hutfnb%$k ! ! "nsta connectpoints"

zone 3 1
n/plot nt.(tf-hutf)/hutf%$k ! ! "nsta connectpoints"
n/plot nt.(tfc-hutfc)/hutfc%$k ! ! "nsta connectpoints"
n/plot nt.(tfb-hutfb)/hutfb%$k ! ! "nsta connectpoints"
zone 2 1
n/plot nt.(tfnw-hutfnw)/hutfnw%$k ! ! "nsta connectpoints"
n/plot nt.(tfnb-hutfnb)/hutfnb%$k ! ! "nsta connectpoints"

######################
zone
n/plot nt.tf%$k ! ! "nsta connectpoints"
n/plot nt.tfnw%$k ! ! "nsta connectpoints red same"
n/plot nt.tfnb%$k ! ! "nsta connectpoints blue same"

zone
n/plot nt.tf/tfnw%$k tfnw>0. ! "nsta connectpoints"
n/plot nt.hutf/hutfnw%$k hutfnw>0. ! "nsta connectpoints same blue"
addline -1000 1 1000 1 "red"

zone
n/plot nt.tf/tfnb%$k tfnb>0. ! "nsta connectpoints"
n/plot nt.hutf/hutfnb%$k hutfnb>0. ! "nsta connectpoints same blue"
addline -1000 1 1000 1 "red"

######################
zone
n/scan nt ! ! -f:cmvchkwhu.data k tf tfnb
*/