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

void compute_xroot(void);

// cmvtstblack [T,nutest(Hz)] [deriv 0/1]

int main(int narg,char *arg[])
{
  //{compute_xroot(); return -41;}
 double T = 2.718, nutest=-1.;
 unsigned short funcderiv = 0;
 double x0;

 if(narg>1) sscanf(arg[1],"%lf,%lf",&T,&nutest);
 if(T<=0.) T = 2.718;
 if(narg>2) sscanf(arg[2],"%hu",&funcderiv);
 PlanckSpectra::SpectraFunc deriv = (funcderiv==0) ? PlanckSpectra::VALUE : PlanckSpectra::DERIV;

 cout<<"Temperature "<<T<<" K , deriv="<<(int)deriv<<endl;

 //--------------------------------- Radiance Frequence
 // frequence planck radiance W/..
 PlanckSpectra pradf(T);
 pradf.SetSpectraFunc(deriv);
 pradf.SetSpectraVar(PlanckSpectra::NU);
 pradf.SetSpectraApprox(PlanckSpectra::PLANCK);
 pradf.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 pradf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phpradf(pradf);
 phpradf.SetSpectraPower(PlanckSpectra::PHOTON);

 // frequence rayleigh radiance W/..
 PlanckSpectra rradf(T);
 rradf.SetSpectraFunc(deriv);
 rradf.SetSpectraVar(PlanckSpectra::NU);
 rradf.SetSpectraApprox(PlanckSpectra::RAYLEIGH);
 rradf.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 rradf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phrradf(rradf);
 phrradf.SetSpectraPower(PlanckSpectra::PHOTON);

 // frequence wien radiance W/..
 PlanckSpectra wradf(T);
 wradf.SetSpectraFunc(deriv);
 wradf.SetSpectraVar(PlanckSpectra::NU);
 wradf.SetSpectraApprox(PlanckSpectra::WIEN);
 wradf.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 wradf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phwradf(wradf);
 phwradf.SetSpectraPower(PlanckSpectra::PHOTON);

 // Check
 if(nutest<=0.) nutest = k_Boltzman_Cst*T/h_Planck_Cst;
 cout<<"Test for t="<<T<<" K,  nu="<<nutest<<" Hz"<<endl
     <<"    planck = "<<pradf(nutest)<<endl
     <<"      wien = "<<wradf(nutest)<<endl
     <<"  rayleigh = "<<rradf(nutest)<<endl;

 //--------------------------------- Densite Frequence
 // frequence planck density W/..
 PlanckSpectra pdensf(T);
 pdensf.SetSpectraFunc(deriv);
 pdensf.SetSpectraVar(PlanckSpectra::NU);
 pdensf.SetSpectraApprox(PlanckSpectra::PLANCK);
 pdensf.SetSpectraUnit(PlanckSpectra::DENSENERG);
 pdensf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phpdensf(pdensf);
 phpdensf.SetSpectraPower(PlanckSpectra::PHOTON);

 // frequence rayleigh density W/..
 PlanckSpectra rdensf(T);
 rdensf.SetSpectraFunc(deriv);
 rdensf.SetSpectraVar(PlanckSpectra::NU);
 rdensf.SetSpectraApprox(PlanckSpectra::RAYLEIGH);
 rdensf.SetSpectraUnit(PlanckSpectra::DENSENERG);
 rdensf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phrdensf(rdensf);
 phrdensf.SetSpectraPower(PlanckSpectra::PHOTON);

 // frequence wien density W/..
 PlanckSpectra wdensf(T);
 wdensf.SetSpectraFunc(deriv);
 wdensf.SetSpectraVar(PlanckSpectra::NU);
 wdensf.SetSpectraApprox(PlanckSpectra::WIEN);
 wdensf.SetSpectraUnit(PlanckSpectra::DENSENERG);
 wdensf.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phwdensf(wdensf);
 phwdensf.SetSpectraPower(PlanckSpectra::PHOTON);

 //--------------------------------- Radiance Longeur d'onde
 // longueur d'onde planck radiance W/..
 PlanckSpectra pradl(T);
 pradl.SetSpectraFunc(deriv);
 pradl.SetSpectraVar(PlanckSpectra::LAMBDA);
 pradl.SetSpectraApprox(PlanckSpectra::PLANCK);
 pradl.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 pradl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phpradl(pradl);
 phpradl.SetSpectraPower(PlanckSpectra::PHOTON);

 // longueur d'onde rayleigh radiance W/..
 PlanckSpectra rradl(T);
 rradl.SetSpectraFunc(deriv);
 rradl.SetSpectraVar(PlanckSpectra::LAMBDA);
 rradl.SetSpectraApprox(PlanckSpectra::RAYLEIGH);
 rradl.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 rradl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phrradl(rradl);
 phrradl.SetSpectraPower(PlanckSpectra::PHOTON);

 // longueur d'onde wien radiance W/..
 PlanckSpectra wradl(T);
 wradl.SetSpectraFunc(deriv);
 wradl.SetSpectraVar(PlanckSpectra::LAMBDA);
 wradl.SetSpectraApprox(PlanckSpectra::WIEN);
 wradl.SetSpectraUnit(PlanckSpectra::ANGSFLUX);
 wradl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phwradl(wradl);
 phwradl.SetSpectraPower(PlanckSpectra::PHOTON);

 //--------------------------------- Densite Longeur d'onde
 // longueur d'onde planck density W/..
 PlanckSpectra pdensl(T);
 pdensl.SetSpectraFunc(deriv);
 pdensl.SetSpectraVar(PlanckSpectra::LAMBDA);
 pdensl.SetSpectraApprox(PlanckSpectra::PLANCK);
 pdensl.SetSpectraUnit(PlanckSpectra::DENSENERG);
 pdensl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phpdensl(pdensl);
 phpdensl.SetSpectraPower(PlanckSpectra::PHOTON);

 // longueur d'onde rayleigh density W/..
 PlanckSpectra rdensl(T);
 rdensl.SetSpectraFunc(deriv);
 rdensl.SetSpectraVar(PlanckSpectra::LAMBDA);
 rdensl.SetSpectraApprox(PlanckSpectra::RAYLEIGH);
 rdensl.SetSpectraUnit(PlanckSpectra::DENSENERG);
 rdensl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phrdensl(rdensl);
 phrdensl.SetSpectraPower(PlanckSpectra::PHOTON);

 // longueur d'onde wien density W/..
 PlanckSpectra wdensl(T);
 wdensl.SetSpectraFunc(deriv);
 wdensl.SetSpectraVar(PlanckSpectra::LAMBDA);
 wdensl.SetSpectraApprox(PlanckSpectra::WIEN);
 wdensl.SetSpectraUnit(PlanckSpectra::DENSENERG);
 wdensl.SetSpectraPower(PlanckSpectra::POWER);
 PlanckSpectra phwdensl(wdensl);
 phwdensl.SetSpectraPower(PlanckSpectra::PHOTON);

 //--------------------------------- Le maximum des spectres de planck
 cout<<endl;
 double eps = 0.001;
 x0 = pradf.WienLaw();
 printf("Planck maximum for radiance energy: f = %e Hz  -> %e W/m^2/sr/Hz\n",x0,pradf(x0));
 printf("                             check:   = %e\n",pradf.FindMaximum(eps));
 x0 = phpradf.WienLaw();
 printf("Planck maximum for radiance photon: f = %e Hz  -> %e ph/s/m^2/sr/Hz\n",x0,phpradf(x0));
 printf("                             check:   = %e\n",phpradf.FindMaximum(eps));

 x0 = pdensf.WienLaw();
 printf("Planck maximum for density energy:  f = %e Hz  -> %e J/m^3/Hz\n",x0,pdensf(x0));
 printf("                             check:   = %e\n",pdensf.FindMaximum(eps));
 x0 = phpdensf.WienLaw();
 printf("Planck maximum for density photon:  f = %e Hz  -> %e ph/m^3/Hz\n",x0,phpdensf(x0));
 printf("                             check:   = %e\n",phpdensf.FindMaximum(eps));

 x0 = pradl.WienLaw();
 printf("Planck maximum for radiance energy: l = %e m  -> %e W/m^2/sr/m\n",x0,pradl(x0));
 printf("                             check:   = %e\n",pradl.FindMaximum(eps));
 x0 = phpradl.WienLaw();
 printf("Planck maximum for radiance photon: l = %e m  -> %e ph/s/m^2/sr/m\n",x0,phpradl(x0));
 printf("                             check:   = %e\n",phpradl.FindMaximum(eps));

 x0 = pdensl.WienLaw();
 printf("Planck maximum for density energy:  l = %e m  -> %e J/m^3/m\n",x0,pdensl(x0));
 printf("                             check:   = %e\n",pdensl.FindMaximum(eps));
 x0 = phpdensl.WienLaw();
 printf("Planck maximum for density photon:  l = %e m  -> %e ph/m^3/m\n",x0,phpdensl(x0));
 printf("                             check:   = %e\n",phpdensl.FindMaximum(eps));

 //--------------------------------- La loi de Wien
 cout<<endl;
 x0 = wradf.WienLaw();
 printf("Wien law for radiance energy: f = %e Hz  -> %e W/m^2/sr/Hz\n",x0,wradf(x0));
 printf("                       check:   = %e\n",wradf.FindMaximum(eps));
 x0 = phwradf.WienLaw();
 printf("Wien law for radiance photon: f = %e Hz  -> %e ph/s/m^2/sr/Hz\n",x0,phwradf(x0));
 printf("                       check:   = %e\n",phwradf.FindMaximum(eps));

 x0 = wdensf.WienLaw();
 printf("Wien law for density energy:  f = %e Hz  -> %e J/m^3/Hz\n",x0,wdensf(x0));
 printf("                      check:    = %e\n",wdensf.FindMaximum(eps));
 x0 = phwdensf.WienLaw();
 printf("Wien law for density photon:  f = %e Hz  -> %e  ph/m^3/Hz\n",x0,phwdensf(x0));
 printf("                      check:    = %e\n",phwdensf.FindMaximum(eps));

 x0 = wradl.WienLaw();
 printf("Wien law for radiance energy: l = %e m  -> %e W/m^2/sr/m\n",x0,wradl(x0));
 printf("                       check:   = %e\n",wradl.FindMaximum(eps));
 x0 = phwradl.WienLaw();
 printf("Wien law for radiance photon: l = %e m  -> %e ph/s/m^2/sr/m\n",x0,phwradl(x0));
 printf("                       check:   = %e\n",phwradl.FindMaximum(eps));

 x0 = wdensl.WienLaw();
 printf("Wien law for density energy:  l = %e m  -> %e J/m^3/m\n",x0,wdensl(x0));
 printf("                      check:    = %e\n",wdensl.FindMaximum(eps));
 x0 = phwdensl.WienLaw();
 printf("Wien law for density photon:  l = %e m  -> %e ph/m^3/m\n",x0,phwdensl(x0));
 printf("                      check:    = %e\n",phwdensl.FindMaximum(eps));

 //--------------------------------- Densite totale d'energie
 cout<<endl;
 cout<<"Densite totale d'energie: "<<pdensf.PlanckEnergie()<<" J/m^3"<<endl;
 if(! deriv) {
   double x1,x2,perc,dxinc;
   x0 = pdensf.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
   cout<<"   check by integration: "<<IntegrateFunc(pdensf,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
   cout<<"   check by integration: "<<IntegrateFuncLog(pdensf,x1,x2,perc,dxinc)<<" J/m^3"<<endl;

   x0 = pdensl.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
   cout<<"   check by integration: "<<IntegrateFunc(pdensl,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
   cout<<"   check by integration: "<<IntegrateFuncLog(pdensl,x1,x2,perc,dxinc)<<" J/m^3"<<endl;
 }

 //--------------------------------- Densite totale de photons
 cout<<endl;
 cout<<"Densite totale de photons: "<<pdensf.PlanckPhoton()<<" ph/m^3"<<endl;
 if(! deriv) {
   double x1,x2,perc,dxinc;
   x0 = phpdensf.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
   cout<<"   check by integration: "<<IntegrateFunc(phpdensf,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
   cout<<"   check by integration: "<<IntegrateFuncLog(phpdensf,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;

   x0 = phpdensl.WienLaw(); x1=x0/100.; x2=x0*100.; perc=0.1; dxinc=x0/1000.;
   cout<<"   check by integration: "<<IntegrateFunc(phpdensl,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
   x1=log10(x0/100.); x2=log10(x0*100.); perc=0.1; dxinc=(x2-x1)/100.;
   cout<<"   check by integration: "<<IntegrateFuncLog(phpdensl,x1,x2,perc,dxinc)<<" ph/m^3"<<endl;
 }

 //--------------------------------- NTuple
 cout<<endl;
 int npt = 1000;
 double frac = 10000.;
 x0 = wradf.WienLaw();
 double xmin = x0/frac, xmax = x0*frac;
 cout<<"xmin="<<xmin<<" Hz,  xmax="<<xmax<<" Hz"<<endl;
 double lnx1 = log10(xmin), lnx2 = log10(xmax), dlnx = (lnx2-lnx1)/npt;
 cout<<"lnx1="<<lnx1<<",  lnx2="<<lnx2<<", dlnx="<<dlnx<<endl;

 const int n = 14;
 const char *vname[n] = {"f","l","prf","prl","pdf","pdl","rrf","rrl","rdf","rdl","wrf","wrl","wdf","wdl"};
 NTuple nt(n,vname);
 double xnt[n];
 for(double lx=lnx1;lx<lnx2+dlnx/2.;lx+=dlnx) {
   double f = pow(10.,lx);
   double l = wradf.F2L(f);
   xnt[0] = f;
   xnt[1] = l;
   xnt[2] = pradf(f);
   xnt[3] = pradl(l);
   xnt[4] = pdensf(f);
   xnt[5] = pdensl(l);
   xnt[6] = rradf(f);
   xnt[7] = rradl(l);
   xnt[8] = rdensf(f);
   xnt[9] = rdensl(l);
   xnt[10] = wradf(f);
   xnt[11] = wradl(l);
   xnt[12] = wdensf(f);
   xnt[13] = wdensl(l);
   nt.Fill(xnt);
 }

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

 return 0;
}

//---------------------------------------
void compute_xroot(void)
// Pour calculer les x qui donnent le maximum de f(x) pour:
// f(x) = x^n * exp(-x)
// ou
// f(x) = x^n / (exp(x)-1)
// ou
// f(x) = x^n * exp(x) / (exp(x)-1)^2
{
  double n = 3.;
  cout<<"n="<<n<<endl;
  int npt=10;
  double xmin=0.1,xmax=10., dx=(xmax-xmin)/npt;
  for(int itry=0;itry<100;itry++) {
    double x0=xmin, v0=-1e10;
    for(double x=xmin;x<xmax+dx/2.; x+=dx) {
      double f = pow(x,n) * exp(-x);  // Wien ou dWein/dT
      //double f = pow(x,n) / (exp(x)-1);  // Planck
      //double f = pow(x,n)  * exp(x) / ((exp(x)-1)*(exp(x)-1));  // dPlanck/dT
      if(f>v0) {v0=f; x0=x;}
    }
    printf("x0 = %.18f , v0=%.18f , dx=%e\n",x0,v0,dx);
    xmin = x0-dx; xmax = x0+dx; dx=(xmax-xmin)/npt;
    if(dx<1e-15) break;
  }
}


/*
openppf cmvtstblack.ppf

# Energie
set fac 1.
# Photon
set fac (6.6260693e-34*f)

n/plot nt.log10(f)%log10(l)

# radiance frequence
n/plot nt.prf/${fac}%log10(f) ! ! "nsta connectpoints"
n/plot nt.rrf/${fac}%log10(f) ! ! "nsta connectpoints red same"
n/plot nt.wrf/${fac}%log10(f) ! ! "nsta connectpoints blue same"

# radiance longeur d'onde
n/plot nt.prl/${fac}%log10(l) ! ! "nsta connectpoints"
n/plot nt.rrl/${fac}%log10(l) ! ! "nsta connectpoints red same"
n/plot nt.wrl/${fac}%log10(l) ! ! "nsta connectpoints blue same"

# densite frequence
n/plot nt.pdf/${fac}%log10(f) ! ! "nsta connectpoints"
n/plot nt.rdf/${fac}%log10(f) ! ! "nsta connectpoints red same"
n/plot nt.wdf/${fac}%log10(f) ! ! "nsta connectpoints blue same"

# densite longeur d'onde
n/plot nt.pdl/${fac}%log10(l) ! ! "nsta connectpoints"
n/plot nt.rdl/${fac}%log10(l) ! ! "nsta connectpoints red same"
n/plot nt.wdl/${fac}%log10(l) ! ! "nsta connectpoints blue same"

# Check:  p(l)dl = p(f)df avec dl = c/f^2 df = l^2/c df
n/plot nt.prf%log10(f) ! ! "nsta connectpoints"
n/plot nt.prl*2.998e8/(f*f)%log10(f) ! ! "nsta connectpoints red same"

 */