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

// cmvtvarspec [zref,R]
void usage(void);
void usage(void)
{
  cout<<"cmvtvarspec -z zref -R R -n as,ns -k kmin,kmax -i eps,npt,fracmax"<<endl;
}


int main(int narg,char *arg[])
{
 double ob0 = 0.0444356;
 // -- WMAP
 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
 // -- ouvert matter only
 //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
 // -- plat matter only
 //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;

 double ns = 1., as = 1.;

 double zref = 0.5;
 double R=8./h100;

 double kmin=1e-5,kmax=1000.;
 int npt = 10000;
 double eps=1.e-3, fracmax=1e-4;

 char c;
 while((c = getopt(narg,arg,"hz:R:n:k:i:")) != -1) {
  switch (c) {
  case 'z' :
    sscanf(optarg,"%lf",&zref);
    break;
  case 'R' :
    sscanf(optarg,"%lf",&R);
    if(R<0.) R = -R/h100;
    break;
  case 'n' :
    sscanf(optarg,"%lf,%lf",&as,&ns);
    break;
  case 'k' :
    sscanf(optarg,"%lf,%lf",&kmin,&kmax);
    break;
  case 'i' :
    sscanf(optarg,"%lf,%d,%lf",&eps,&npt,&fracmax);
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }

 double Rg=R/sqrt(5.);
 double lkmin=log10(kmin), lkmax=log10(kmax), dlk=(lkmax-lkmin)/npt;

 cout<<"zref="<<zref<<endl;
 cout<<"as="<<as<<" ns="<<ns<<endl;
 cout<<"R="<<R<<", Rg="<<Rg<<endl;
 cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<")"
     <<", dlk="<<dlk<<" npt="<<npt<<endl;
 cout<<"eps="<<eps<<"  fracmax="<<fracmax<<endl;

 //-----------------------------------------------------------------
 InitialSpectrum pkini(ns,as);

 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
 //tf.SetNoOscEnv(2);

 GrowthFactor d1(om0,ol0);

 PkSpectrum0 pk0(pkini,tf);

 PkSpectrumZ pkz(pk0,d1,zref);

 //-----------------------------------------------------------------
 double k1,k2;
 int rc = 0;
 //----
 cout<<endl<<"Filtrage top_hat"<<endl;
 VarianceSpectrum varpk_th(pkz,0);
 double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
 double pkmax_th = varpk_th(kfind_th);
 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
 k1=kmin, k2=kmax;
 rc = varpk_th.FindLimits(R,pkmax_th*fracmax,k1,k2,eps);
 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;

 varpk_th.SetInteg(0.01,dlk,-1.,4);
 cout<<"varpk_th="<<varpk_th.Variance(R,k1,k2)<<endl;

 //----
 cout<<endl<<"Filtrage gaussien"<<endl;
 VarianceSpectrum varpk_ga(pkz,1);
 double kfind_ga = varpk_ga.FindMaximum(Rg,kmin,kmax,eps);
 double pkmax_ga = varpk_ga(kfind_ga);
 cout<<"kfind_ga = "<<kfind_ga<<" -> "<<log10(kfind_ga)<<", integrand="<<pkmax_ga<<endl;
 k1=kmin, k2=kmax;
 rc = varpk_ga.FindLimits(Rg,pkmax_ga*fracmax,k1,k2,eps);
 cout<<"limit_ga: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;

 varpk_ga.SetInteg(0.01,dlk,-1.,4);
 cout<<"varpk_ga="<<varpk_ga.Variance(Rg,k1,k2)<<endl;
 
 //----
  cout<<endl<<"Filtrage 1 (integrale du sppectre)"<<endl;
 VarianceSpectrum varpk_int(pkz,2);
 double kfind_int = varpk_int.FindMaximum(Rg,kmin,kmax,eps);
 double pkmax_int = varpk_int(kfind_int);
 cout<<"kfind_int = "<<kfind_int<<" -> "<<log10(kfind_int)<<", integrand="<<pkmax_int<<endl;
 k1=kmin, k2=kmax;
 rc = varpk_int.FindLimits(Rg,pkmax_int*fracmax,k1,k2,eps);
 cout<<"limit_int: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "<<k2<<" ("<<log10(k2)<<")"<<endl;

 varpk_int.SetInteg(0.01,dlk,-1.,4);
 cout<<"varpk_int="<<varpk_int.Variance(Rg,k1,k2)<<endl;
 
//-----------------------------------------------------------------
 const int n = 2;
 char *vname[n] = {"k","pk"};
 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] = pkz(k);
   nt.Fill(xnt);
 }

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

 return 0;
}

/*
openppf cmvtvarspec.ppf

#### Spectre a z
n/plot nt.pk%log10(k) ! ! "nsta crossmarker3"

# L'approximation pour j1(x) (top-hat) et la gaussienne equivalente
set l1 1e-6
set l2 10
set j1sx pow(3*(sin(x)-x*cos(x))/(x*x*x),2.)
set j1sxa3 (1.-x*x/5.)
set j1sxa5 (1.-x*x/5.*(1.-3.*x*x/35.))
set j1sxa7 (1.-x*x/5.*(1.-3.*x*x/35.*(1.-4.*x*x/81.)))
func $j1sx    $l1 $l2 1000
func $j1sxa3  $l1 $l2 1000 "red same"
func $j1sxa5  $l1 $l2 1000 "blue same"
func $j1sxa7  $l1 $l2 1000 "orange same"
set gauss  pow(exp(-(x/sqrt(5.))*(x/sqrt(5.))/2.),2.)
func $gauss  $l1 $l2 1000 "green same"

func $j1sx-$j1sxa3  $l1 $l2 1000 "red"
func $j1sx-$j1sxa5  $l1 $l2 1000 "blue same"
func $j1sx-$j1sxa7  $l1 $l2 1000 "orange same"
func $j1sx-$gauss  $l1 $l2 1000 "green same"
addline $l1 0 $l2 0

# la fonction "Pk" a integrer
zone
set D (k*k*pk/(2.*M_PI*M_PI))
n/plot nt.$D%log10(k) ! ! "nsta crossmarker3"

# La fonction a integrer en dk lineaire (top-hat et gaussien)
set R 10
set RG ($R/sqrt(5.))
set x (k*$R)
set j1sx pow(3*(sin($x)-$x*cos($x))/($x*$x*$x),2.)
set y (k*$RG)
set gauss  pow(exp(-$y*$y/2.),2.)

n/plot nt.$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
n/plot nt.$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"

# La fonction a integrer en d(log(x)) (top-hat et gaussien)
# c'est ce que fait Variance()
n/plot nt.k*$D*$j1sx%log10(k) ! ! "nsta crossmarker3"
n/plot nt.k*$D*$gauss%log10(k) ! ! "nsta crossmarker3 same red"
n/plot nt.k*$D%log10(k) ! ! "nsta crossmarker3 same blue"

*/