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

void usage(void);
void usage(void) {cout<<"cmvtgrowth z1,z2,nz [Omatter0,Lambda0]"<<endl;}

int main(int narg,char *arg[])
{

 double z1=0., z2=10.;
 int nz = 100;;
 if(narg>1) sscanf(arg[1],"%lf,%lf,%d",&z1,&z2,&nz);
 if(nz<=0) nz = 100;
 double dz = (z2-z1)/nz;
 cout<<"z1="<<z1<<"  z2="<<z2<<"  nz="<<nz<<"  dz="<<dz<<endl;

 double om0=0.267804, ol0=0.73;
 if(narg>2) sscanf(arg[2],"%lf,%lf",&om0,&ol0);
 cout<<"Om0="<<om0<<"  Ol0="<<ol0<<endl;

 GrowthEisenstein growth(om0,ol0);
 cout<<"D1(z=0) = "<<growth(0.)<<endl;

 const int n = 3;
 const char *vname[n] = {"z","d1","d1dz"};
 NTuple nt(n,vname);
 double xnt[n];
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   xnt[0]  = z;
   xnt[1]  = growth(z);
   xnt[2]  = growth.DsDz(z,dz/2.);
   nt.Fill(xnt);
 }

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

 return 0;
}


/*
openppf cmvtgrowth.ppf

set cut 1
set cut z<5

# --- growth
n/plot nt.d1%z $cut ! "nsta connectpoints"
n/plot nt.1./(1.+z)%z $cut ! "nsta connectpoints red same"

# --- growth'/growth
n/plot nt.d1dz/d1%z $cut ! "nsta connectpoints"
n/plot nt.-1./(1.+z)%z $cut ! "nsta connectpoints red same"

# --- d(growth)/dz
zmin = 0.
zmax = 10.
set npt ${dv.size}
dd = ($zmax-$zmin)/(${npt}-1.)
exptovec dv nt d1
c++exec for(int i=1;i<dv.Size(); i++) dv(i-1) = dv(i)-dv(i-1);

n/plot nt.-1./pow(1.+z,2.)%z $cut ! "nsta connectpoints red"
n/plot nt.d1dz%z $cut ! "nsta connectpoints same"
n/plot dv.val/${dd}%${zmin}+(n+1)*$dd ! ! "nsta plusmarker5 green same"

*/