#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
<<"cmvtransf [options]"<<endl
<<"  -k k1,k2,npt : k range in Mpc^-1"<<endl
<<"  -U h100,Om0,Ob0,tcmb : cosmology"<<endl
<<"  -F filename : read also CMBFast file"<<endl
<<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;
 string fcmbfast = "";

 char c;
  while((c = getopt(narg,arg,"hk:U:F:")) != -1) {
  switch (c) {
  case 'k' :
    sscanf(optarg,"%d,%lf,%lf",&npt,&k1,&k2);
    break;
  case 'U' :
    sscanf(optarg,"%lf,%lf,%lf,%lf",&h100,&Om0,&Ob0,&tcmb);
    break;
  case 'F' :
    fcmbfast = optarg;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }

 if(k1<=0.) {cout<<"k1 must be >0"<<endl; return -1;}
 if(npt<=0) npt = 1000;
 if(k2<=k1) {k1=1.e-4; k2=1.e+2;}
 cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
 cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
 cout<<"fcmbfast="<<fcmbfast<<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

 TransfertTabulate tfcmbfast;
 bool cmbfastOK = false;
 if(fcmbfast.size()>0) {
   cout<<endl<<"TransfertTabulate for CMBfast"<<endl;
   tfcmbfast.ReadCMBFast(fcmbfast,h100,Om0-Ob0,Ob0);
   tfcmbfast.SetInterpTyp(1);
   if(tfcmbfast.NPoints()>0) cmbfastOK = true;
 }

 const int n = 6;
 const char *vname[n] = {"k","t","tnosc1","tnosc2","tnob","tcmbf"};
 NTuple nt(n,vname);
 double xnt[n];
 for(int i=0;i<n;i++) xnt[i]=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);
   xnt[0]  = k;
   xnt[1]  = tf(k);
   xnt[2]  = tfnosc1(k);
   xnt[3]  = tfnosc2(k);
   xnt[4]  = tfnob(k);
   if(cmbfastOK) xnt[5]  = tfcmbfast(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)

# Eisenstein
n/plot nt.t%$k ! ! "nsta connectpoints"
n/plot nt.tnob%$k ! ! "nsta red same connectpoints"
n/plot nt.tnosc2%$k ! ! "nsta blue same connectpoints"
n/plot nt.tnosc1%$k ! ! "nsta green same connectpoints"

n/plot nt.t/tnob%$k   tnob>0   ! "nsta red connectpoints"
n/plot nt.t/tnosc2%$k tnosc2>0 ! "nsta blue same connectpoints"
n/plot nt.t/tnosc1%$k tnosc1>0 ! "nsta green same connectpoints"

# CMBFast
n/plot nt.t%$k ! ! "nsta connectpoints"
n/plot nt.tcmbf%$k ! ! "nsta red same connectpoints"

n/plot nt.(tcmbf-t)%$k tcmbf>0. ! "nsta connectpoints"
n/plot nt.(tcmbf-t)/tcmbf%$k tcmbf>0. ! "nsta connectpoints"
*/