#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 "luc.h"

void usage(void);
void usage(void)
{
  cout<<"cmvtuniv z1,z2,dz"<<endl;
}

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

 double z1=0., z2=10., dz=1.;
 if(narg>1) sscanf(arg[1],"%lf,%lf,%lf",&z1,&z2,&dz);

 SimpleUniverse univ;
 univ.SetPrtDbgLevel(1);

 univ.Seth(h100);
 univ.SetOmegaMatter(om0);
 univ.SetOmegaLambda(ol0);
 //univ.SetDarkEnergy(ol0,w0);
 univ.SetOmegaRadiation(or0);
 univ.SetOmegaBaryon(ob0);
 univ.SetOmegaPhoton(oph0);

 //univ.SetFlatUniverse_OmegaMatter();
 //univ.SetFlatUniverse_OmegaLambda();

 /*
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   univ.SetEmissionRedShift(z,0.001);
   cout<<endl;
   univ.Print();
 }
 */

 /**/
 double norm = univ.HubbleLengthMpc(); // 1.;
 double norm3 = pow(norm,3.);
 const int n = 11;
 char *vname[n] = {
   "z","hz","om","or","ol","ok","ot",
   "da","dl","dvc","vc0"
 };
 NTuple nt(n,vname);
 double xnt[n];
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   univ.SetEmissionRedShift(z,0.001);
   xnt[0]  = univ.ZE();
   xnt[1]  = univ.HZE();
   xnt[2]  = univ.OmegaMatterZE();
   xnt[3]  = univ.OmegaRadiationZE();
   xnt[4]  = 0.; //univ.OmegaLambdaZE()+univ.OmegaDarkZE();
   xnt[5]  = univ.OmegaCurvZE();
   xnt[6]  = 0.; //univ.OmegaTotalZE();
   xnt[7]  = univ.AngularDiameterDistanceMpc()/norm;
   xnt[8]  = univ.LuminosityDistanceMpc()/norm;
   if(xnt[8]<0.) cout<<xnt[8]<<endl;
   xnt[9]  = 0.; //univ.dVolZE()/norm3;
   xnt[10] = 0.; //univ.Vol4PiZE()/norm3;
   nt.Fill(xnt);
 }
 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvtluc.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
 /**/

 return 0;
}


/*
openppf cmvtluc.ppf

set cut 1
set cut z<100.

zone
n/plot nt.hz%z $cut ! "nsta"

zone
n/plot nt.dl%z $cut ! "nsta"
n/plot nt.da%z $cut ! "nsta same red"

*/