#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 "constcosmo.h"
#include "pkspectrum.h"

void usage(void);
void usage(void) {cout<<"cmvtuniv [k1,k2,npt] [h100,Om0,Ob0,tcmb]"<<endl;}

int main(int narg,char *arg[])
{
 double k1 = 1e-6, k2 = 10.; int_4 npt = 100;
 double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;

 if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&k1,&k2,&npt);
 if(k1<=0.) {cout<<"k1 must be >0"<<endl; return -1;}
 if(npt<=0) npt = 100;
 cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
 if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;

 cout<<"TransfertEisenstein with baryon"<<endl;
 TransfertEisenstein tf(h100,Om0-Ob0,Ob0,tcmb);
 cout<<"KPeak = "<<tf.KPeak()<<"  1/Mpc"<<endl;

 cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.1"<<endl;
 TransfertEisenstein tfnosc1(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.1
 tfnosc1.SetNoOscEnv(1);

 cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.2"<<endl;
 TransfertEisenstein tfnosc2(h100,Om0-Ob0,Ob0,tcmb);  // No oscillation approx.2
 tfnosc2.SetNoOscEnv(2);

 cout<<endl<<"TransfertEisenstein without baryon"<<endl;
 TransfertEisenstein tfnob(h100,Om0-Ob0,Ob0,tcmb,true);  // No baryons

 const int n = 5;
 char *vname[n] = {"k","t","tnosc1","tnosc2","tnob"};
 NTuple nt(n,vname);
 double xnt[n];

 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);
   xnt[0]  = k;
   xnt[1]  = tf(k);
   xnt[2]  = tfnosc1(k);
   xnt[3]  = tfnosc2(k);
   xnt[4]  = tfnob(k);
   nt.Fill(xnt);
 }

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

 return 0;
}

/*
openppf cmvtransf.ppf

set k k
set k log10(k)

n/plot nt.t%$k ! ! "nsta crossmarker3"
n/plot nt.tnob%$k ! ! "nsta red crossmarker3 same"
n/plot nt.tnosc1%$k ! ! "nsta blue circlemarker3 same"
n/plot nt.tnosc2%$k ! ! "nsta green trianglemarker3 same"

n/plot nt.t/tnosc1%$k ! ! "nsta crossmarker3"
n/plot nt.t/tnosc2%$k ! ! "nsta crossmarker3 red same"
*/