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

void usage(void);
void usage(void) {cout<<"cmvtuniv z1,z2,dz [perc,dzinc,dzmax,glorder]"<<endl;}

void tstprint(CosmoCalc& univ,double z1,double z2,double dz);
void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz);
void tstspeed(CosmoCalc& univ,double z1,double z2,double dz);
void tstntuple(CosmoCalc& univ,double z1,double z2,double dz);
void tstinterp(CosmoCalc& univ,double z1,double z2,double dz);

const double deg2rad = M_PI/180.;

int main(int narg,char *arg[])
{
 unsigned short flat = 0;
 // -- 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.; flat = 1;
 // -- plat lambda only
 //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.; flat = 2;

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

 double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
 if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%hu",&perc,&dzinc,&dzmax,&glorder);
 cout<<" perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;

 CosmoCalc univ(flat,true,z2);
 univ.SetInteg(perc,dzinc,dzmax,glorder);
 univ.SetDynParam(h100,om0,or0,ol0,w0);
 univ.PrtInteg();
 univ.Print();

 CosmoCalc univ1(univ);
 univ1.SetInteg(perc,dzinc,dzmax,glorder);
 univ1.SetDynParam(h100,om0,or0,ol0,w0);
 univ1.PrtInteg();
 univ1.Print();

 tstprint(univ,z1,z2,dz);
 //tstprint(univ,univ1,z1,z2,dz);
 //tstspeed(univ,z1,z2,dz);
 //tstntuple(univ,z1,z2,dz);
 //tstinterp(univ,z1,z2,dz);

 return 0;
}

/* ----------------------------------------------------- */
void tstprint(CosmoCalc& univ,double z1,double z2,double dz)
{
 cout<<"\nTSTPRINT()"<<endl;
 double dzz = dz; if(dzz>z2-z1) dzz = 0.01;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   cout<<"--"<<endl;
   univ.Print(z);
   cout<<"Volume comoving in [z,z+"<<dzz<<"] for 4Pi sr: "
       <<univ.Vol4Pi(z,z+dzz)<<" Mpc^3"<<endl;
   double dang = univ.Dang(z);
   double a = deg2rad/3600.;
   cout<<"1\" -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
   a = deg2rad/60.;
   cout<<"1\' -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
   a = deg2rad;
   cout<<"1d -> "<<dang*a<<" Mpc = "<<dang*a*(1.+z)<<" Mpc com"<<endl;
   double dloscom = univ.Dhubble() / univ.E(z);
   cout<<"dz=1 -> "<<dloscom<<" Mpc com"<<endl;

 }
}

/* ----------------------------------------------------- */
void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz)
{
 cout<<"\nTSTPRINT()"<<endl;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   cout<<endl<<"--spline:"<<endl;
   univ1.Print(z);
   cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
       <<univ1.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
   cout<<"--integ:"<<endl;
   univ2.Print(z);
   cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
       <<univ2.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
 }
}

/* ----------------------------------------------------- */
void tstspeed(CosmoCalc& univ,double z1,double z2,double dz)
{
 cout<<"\nTSTSPEED()"<<endl;
 double sum = 0.;
 int_4 n=0;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   sum += univ.Dang(z);
   n++;
 }
 if(n>0) cout<<"... "<<sum/n<<endl;
}

/* ----------------------------------------------------- */
void tstntuple(CosmoCalc& univ,double z1,double z2,double dz)
{
 cout<<"\nTSTNTUPLE()"<<endl;
 double norm =  univ.Dhubble();  // 1.
 double norm3 = pow(norm,3.);
 const int n = 15;
 char *vname[n] = {
   "z","hz","om","or","ol","ok","ot","ob",
   "dtc","dlc","da","dl",
   "dvc","vc0","vcdz"
 };
 NTuple nt(n,vname);
 double xnt[n];
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   xnt[0]  = z;
   xnt[1]  = univ.H(z);
   xnt[2]  = univ.Omatter(z);
   xnt[3]  = univ.Orelat(z);
   xnt[4]  = univ.Olambda(z);
   xnt[5]  = univ.Ocurv(z);
   xnt[6]  = univ.Otot(z);
   xnt[7]  = univ.Obaryon(z);
   xnt[8]  = univ.Dtrcom(z)/norm;
   xnt[9]  = univ.Dloscom(z)/norm;
   xnt[10] = univ.Dang(z)/norm;
   xnt[11] = univ.Dlum(z)/norm;
   xnt[12] = univ.dVol(z)/norm3;
   xnt[13] = univ.Vol4Pi(z)/norm3;
   xnt[14] = univ.Vol4Pi(z,z+dz)/norm3;
   nt.Fill(xnt);
 }
 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvtuniv.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
}


/*
openppf cmvtuniv.ppf

set cut 1
set cut z<100.
set cut z<10.

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

zone
n/plot nt.dl%z $cut ! "nsta connectpoints"
n/plot nt.da%z $cut ! "nsta connectpoints same red"
n/plot nt.dtc%z $cut ! "nsta connectpoints same blue"
n/plot nt.dlc%z $cut ! "nsta connectpoints same green"

zone 2 2
n/plot nt.dvc%z $cut ! "nsta connectpoints"
n/plot nt.vc0%z $cut ! "nsta connectpoints"
n/plot nt.vcdz%z $cut ! "nsta connectpoints"

zone
n/plot nt.ot%z $cut ! "nsta connectpoints"
n/plot nt.om%z $cut ! "nsta connectpoints same blue"
n/plot nt.or%z $cut ! "nsta connectpoints same red"
n/plot nt.ol%z $cut ! "nsta connectpoints same green"
n/plot nt.ok%z $cut ! "nsta connectpoints same orange"

 */

/* ----------------------------------------------------- */
void tstinterp(CosmoCalc& univ,double z1,double z2,double dz)
{
 cout<<"\nTSTINTERP()"<<endl;

 // On repmplit les donnes completes
 vector<double> Z,D;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   Z.push_back(z);
   D.push_back(univ.Dloscom(z));
 }
 // On remplit un sous tableau
 int ninc = 5;
 vector<double> Dinterp;
 double zmin = Z[0], zmax = Z[0];
 for(unsigned int i=0;i<Z.size();i+=ninc) {
   Dinterp.push_back(D[i]);
   zmax = Z[i];
 }
 InterpFunc interp(zmin,zmax,Dinterp);
 unsigned short ok;

 const int n = 8;
 char *vname[n] = {"z","d","di","dl","dp","zi","zl","zp"};
 NTuple nt(n,vname);
 double xnt[n];
   
 for(unsigned int i=0;i<Z.size();i++) {
   if(Z[i]>zmax) break;
   xnt[0] = Z[i];
   xnt[1] = D[i];
   xnt[2] = interp(Z[i]);
   xnt[3] = interp.Linear(Z[i],ok);
   xnt[4] = interp.Parab(Z[i],ok);
   xnt[5] = xnt[6] = xnt[7] = 0.;
   nt.Fill(xnt);
 }
 cout<<">>>> Ecriture"<<endl;
 string tag = "cmvtuniv.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
}

/*
openppf cmvtuniv.ppf

n/plot nt.d%z ! ! "nsta connectpoints"
n/plot nt.di%z ! ! "nsta connectpoints same green"
n/plot nt.dl%z ! ! "nsta connectpoints same red"
n/plot nt.dp%z ! ! "nsta connectpoints same blue"

n/plot nt.di-d%z ! ! "nsta connectpoints green"
n/plot nt.dl-d%z ! ! "nsta connectpoints same red"
n/plot nt.dp-d%z ! ! "nsta connectpoints same blue"

n/plot nt.(di-d)/d%z d>0 ! "nsta connectpoints green"
n/plot nt.(dl-d)/d%z d>0 ! "nsta connectpoints same red"
n/plot nt.(dp-d)/d%z d>0 ! "nsta connectpoints same blue"

*/