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

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.;

void usage(void);
void usage(void) {
 cout
 <<"cmvtuniv [options] z1,z2"<<endl
 <<"  -D dz"<<endl
 <<"  -U h100,om0,or0,ol0,w0,wa"<<endl
 <<"  -F flat(0,1,2)"<<endl
 <<"  -I perc,dzinc,dzmax,glorder"<<endl
 <<" --- What has to be done:"<<endl
 <<"  -W 1/2/3/4/5"<<endl
 <<"     1 : print info from z1 to z2 by step dz"<<endl
 <<"     2 : compare with/without spline from z1 to z2 by step dz"<<endl
 <<"     3 : test speed from z1 to z2 by step dz"<<endl
 <<"     4 : fill ntuple with infos from z1 to z2 by step dz"<<endl
 <<"     5 : test interpolation loscom from z1 to z2 by step dz"<<endl
 <<endl;
}

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.,wa=0.;
 // -- ouvert matter only
 //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1,wa=0..;
 // -- plat matter only
 //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.,wa=0.; flat = 1;
 // -- plat lambda only
 //double h100=0.71, om0=0., or0=0., ol0=1.,w0=-1.,wa=0.; flat = 2;

 // -- parametre d'integration
 double perc=0.01,dzinc=-1.,dzmax=-1.;
 unsigned short glorder=4;
 // -- redshift
 double z1=0., z2=-1., dz=-1.;
 // -- Work to be done
 int todo = 1;

 InitTim();

 // Read arguments
 char c;
 while ((c = getopt(narg, arg, "hD:U:F:I:W:")) != -1) {
  switch (c) {
  case 'D' :
    sscanf(optarg,"%lf",&dz);
    break;
  case 'U' :
    sscanf(optarg,"%lf,%lf,%lf,%lf,%lf,%lf",&h100,&om0,&or0,&ol0,&w0,&wa);
    break;
  case 'F' :
    sscanf(optarg,"%hu",&flat);
    break;
  case 'I' :
    sscanf(optarg,"%lf,%lf,%lf,%hu",&perc,&dzinc,&dzmax,&glorder);
    break;
  case 'W' :
    sscanf(optarg,"%d",&todo);
    break;
  case 'h' :
  default:
    usage();
    return -1;
  }
}
 if(optind<narg) sscanf(arg[optind],"%lf,%lf",&z1,&z2);

 if(z1<0.) z1 = 0.;
 if(z2<=z1) z2 = z1+1.;
 if(dz<=0.) dz = 10.*(z2-z1);
 cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
 cout<<"perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
 cout<<"h100="<<h100<<" om0="<<om0<<" or0="<<or0<<" ol0="<<ol0<<" w0="<<w0<<" wa="<<wa<<" flat="<<flat<<endl;
 cout<<"todo ="<<todo<<endl;

 CosmoCalc univ1(flat,true,z2);
 univ1.SetInteg(perc,dzinc,dzmax,glorder);
 univ1.SetDynParam(h100,om0,or0,ol0,w0,wa);
 univ1.PrtInteg();
 univ1.Print();
 cout<<endl;

 if(todo==1)  tstprint(univ1,z1,z2,dz);
 else if(todo==2) {
   CosmoCalc univ2(flat,false,z2);
   univ2.SetInteg(perc,dzinc,dzmax,glorder);
   univ2.SetDynParam(h100,om0,or0,ol0,w0,wa);
   univ2.PrtInteg();
   univ2.Print();
   tstprint(univ1,univ2,z1,z2,dz);
 }
 else if(todo==3) tstspeed(univ1,z1,z2,dz);
 else if(todo==4) tstntuple(univ1,z1,z2,dz);
 else if(todo==5) tstinterp(univ1,z1,z2,dz);
 else cout<<"Unknown job to be done !!!!!  todo="<<todo<<endl;

 return 0;
}

/* ----------------------------------------------------- */
void tstprint(CosmoCalc& univ,double z1,double z2,double dz)
{
 if(dz>z2-z1) dz = z2-z1;
 double fz1 = 1000.*Fr_HyperFin_Par/(1.+z1);
 double fz1pdz = 1000.*Fr_HyperFin_Par/(1.+z1+dz);
 double fz2 = 1000.*Fr_HyperFin_Par/(1.+z2);
 double fz2pdz = 1000.*Fr_HyperFin_Par/(1.+z2+dz);
 cout<<"\nTSTPRINT(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
 cout<<"  "<<fz1<<" -> "<<fz2<<" MHz,    df: "<<fz1-fz1pdz<<" , "<<fz2-fz2pdz<<" MHz"<<endl;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   cout<<"\n-- z = "<<z<<"      ("<<1000.*Fr_HyperFin_Par/(1.+z)<<" MHz)"<<endl;
   univ.Print(z);
   cout<<"Volume comoving in [z,z+"<<dz<<"] for 4Pi sr: "
       <<univ.Vol4Pi(z,z+dz)<<" 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<<"dLOScom/dz -> "<<dloscom<<" Mpc com"<<endl;
 }

}

/* ----------------------------------------------------- */
void tstprint(CosmoCalc& univ1,CosmoCalc& univ2,double z1,double z2,double dz)
{
 if(dz>z2-z1) dz = z2-z1;
 cout<<"\nTSTPRINT(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   cout<<endl<<"--spline: z = "<<z<<endl;
   univ1.Print(z);
   cout<<"Volume comoving in [z,z+dz] for 4Pi sr: "
       <<univ1.Vol4Pi(z,z+dz)<<" Mpc^3"<<endl;
   cout<<"--integ: z = "<<z<<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)
{
 if(dz>z2-z1) dz = z2-z1;
 cout<<"\nTSTSPEED(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
 double sum = 0.;
 int_4 n=0;
 PrtTim("-- Beginning");
 for(double z=z1;z<z2+dz/2.;z+=dz) {
   sum += univ.Dang(z);
   n++;
 }
 if(n>0) cout<<"n="<<n<<" ...dum="<<sum/n<<endl;
 PrtTim("-- Ending");
}

/* ----------------------------------------------------- */
void tstntuple(CosmoCalc& univ,double z1,double z2,double dz)
{
 if(dz>z2-z1) dz = z2-z1;
 cout<<"\nTSTNTUPLE(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
 double norm =  univ.Dhubble();  // 1.
 double norm3 = pow(norm,3.);
 const int n = 15;
 const 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;
 nt.Info()["flat"] = univ.Flat();
 nt.Info()["h100"] = univ.h100();
 nt.Info()["Olambda0"] = univ.Olambda0();
 nt.Info()["W0"] = univ.W0();
 nt.Info()["Wa"] = univ.Wa();
 nt.Info()["Omatter0"] = univ.Omatter0();
 nt.Info()["Obaryon0"] = univ.Obaryon0();
 nt.Info()["Orelat0"] = univ.Orelat0();
 nt.Info()["Ocurv0"] = univ.Ocurv0();
 nt.Info()["Otot0"] = univ.Otot0();
 nt.Info()["ZMax"] = univ.ZMax();

 string tag = "cmvtuniv_nt.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
}


/*
openppf cmvtuniv_nt.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)
{
 if(dz>z2-z1) dz = z2-z1;
 cout<<"\nTSTINTERP(): z1="<<z1<<" z2="<<z2<<" dz="<<dz<<endl;
 // On remplit 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;
 const 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_int.ppf";
 POutPersist pos(tag);
 tag = "nt"; pos.PutObject(nt,tag);
}

/*
openppf cmvtuniv_int.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"

*/