#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "pexceptions.h"

#include "constcosmo.h"
#include "planckspectra.h"

namespace SOPHYA {

/*
----------------------------------------------------------------
- UNITES:
     Temperature en Kelvin  "K"
     Frequence en GigaHertz  "Hz"
     Longueur d'onde en Metre  "m"
----------------------------------------------------------------
- Frequence / Longeur d'onde
 l = c/f ; f = c/l
 dl = c/f^2 df = l^2/c df
 df = c/l^2 dl = f^2/c dl
 P[f] df = P[l] dl = P[f] c/l^2 dl = P[l] c/f^2 df
----------------------------------------------------------------
- Lois de Planck Wien et Rayleigh
     Planck   H*Nu/(K*T) = H*C/(Lambda*K*T) quelconque
     Wien     H*Nu/(K*T) = H*C/(Lambda*K*T) >> 1
     Rayleigh H*Nu/(K*T) = H*C/(Lambda*K*T) << 1
----------------------------------------------------------------
- Les termes en facteur pour les lois de Planck Wien et Rayleigh
        (la partie contenant l'exponentielle)
 ainsi que pour les derivees des ces termes par rapport a la temperature
 Planck:   1/[exp(H*f/(K*t))-1] = 1/[exp(HC/(l*K*t))-1]
     der=  H*f/(K*t) /t * exp(HC/(l*K*t))/[exp(HC/(l*K*t))-1]^2
 Wien:     1/exp(H*f/(K*t))
     der=  H*f/(K*t) /t /exp(H*f/(K*t))
 Rayleigh: 1/(H*f/(K*t))
     der=  K/(H*f)
 Pour x->0: compute exp(x)=1+x+x^2/2 for 1/(exp(x)-1)
            soit 1/(exp(x)-1)  -> 1/(x+x^2/2)
  car exp(x)-1 vaut zero plus vite que son DL a cause de la soustraction
----------------------------------------------------------------
 - RADIANCE (Brillance) en fonction de nu ou lambda
 nombre de watts (ou photons par seconde)
 radies par unite de surface par steradian et par hertz (ou metre)
    2*f^2/c^2   * 1/(exp()-1)    ph/s/m^2/Sr/Hz
    2*h*f^3/c^2 * 1/(exp()-1)    W/m^2/Sr/Hz
    2*c/l^4     * 1/(exp()-1)    ph/s/m^2/Sr/m
    2*h*c^2/l^5 * 1/(exp()-1)    W/m^2/Sr/m
----------------------------------------------------------------
 - EMITTANCE en fonction de nu ou lambda
 nombre de watts (ou photons par seconde)
 radies dans un hemisphere (2*PI sr)
 par unite de surface et par hertz (ou metre)
 ..quand on integre sur 2Pi sr il faut
 ..multiplier par cos(theta)  -> facteur Pi uniquement
    EMITTANCE = Pi * RADIANCE
    2*Pi*f^2/c^2   * 1/(exp()-1)    ph/s/m^2/Hz
    2*Pi*h*f^3/c^2 * 1/(exp()-1)    W/m^2/Hz
    2*Pi*c/l^4     * 1/(exp()-1)    ph/s/m^2/m
    2*Pi*h*c^2/l^5 * 1/(exp()-1)    W/m^2/m
----------------------------------------------------------------
 - DENSITE en fonction de nu ou lambda
 nombre de joules (ou photons) par m^3  et par hertz (ou metre)
   DENSITE = 4*Pi/c * RADIANCE
   8*Pi*f^2/c^3    * 1/(exp()-1)    ph/m^3/Hz
   8*Pi*h*f^3/c^3  * 1/(exp()-1)    J/m^3/Hz
   8*Pi/l^4        * 1/(exp()-1)    ph/m^3/m
   8*Pi*h*c/l^5    * 1/(exp()-1)    J/m^3/m
----------------------------------------------------------------
 - DENSITE d'energie totale  (loi de Stefan-Boltzman)
   4*sigma*T^4/C   J/m^3
   Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
----------------------------------------------------------------
*/

/*
POUR TESTER LES LIMITES DE L'EXPONENTIELLE:
1-/ copier/coller les lignes suivantes dans toto.icc
double x=10.;
for(int i=0;i<35;i++) {
  double ex = exp(x);
  double exm1a = x*(1.+x/2.*(1.+x/3.)),;
  printf("x=%g exp(x)=%e exp(x)-1=%e dl=%e (d=%e)\n",x,ex,ex-1.,exm1a,ex-1.-exm1a);
  x /= 10.;
}
2-/ executer sous spiapp: > c++execfrf toto.icc
*/


#define MAXARGEXP 100.    /* maxi = 709.7827128933840 */
#define MINARGEXP 1.e-10  /* Replace exp(x) by DL for 1/(exp(x)-1) when x->0 */

//************************************************************************
// vitesse de la lumiere en m/s
double CinM_ = SpeedOfLight_Cst*1e+3;
// h/k pour une frequence en Hz
double HsK_ = h_Planck_Cst/k_Boltzman_Cst;
// hc/k pour une longeur d'onde en metre
double HCsK_ = h_Planck_Cst*CinM_/k_Boltzman_Cst;

//************************************************************************
PlanckSpectra::PlanckSpectra(double T)
  : T_(T) , spectraapprox_(PLANCK) , spectrafunc_(VALUE) , spectravar_(NU)
  , spectraunit_(ANGSFLUX) , spectrapower_(POWER)
{
}

PlanckSpectra::PlanckSpectra(PlanckSpectra& s)
  : T_(s.T_) , spectraapprox_(s.spectraapprox_) , spectrafunc_(s.spectrafunc_) , spectravar_(s.spectravar_)
  , spectraunit_(s.spectraunit_), spectrapower_(s.spectrapower_)
{
}

PlanckSpectra::~PlanckSpectra(void)
{
}

//************************************************************************
void PlanckSpectra::SetSpectraApprox(SpectraApprox spectraapprox)
// spectraapprox = PLANCK   : Planck spectra
//               = RAYLEIGH : Rayleigh spectra
//               = WIEN     : Wien spectra
{
  spectraapprox_ = spectraapprox;
}

void PlanckSpectra::SetSpectraFunc(SpectraFunc spectrafunc)
// spectrafunc = VALUE : spectra
//             = DERIV : derivative of spectra relative to temperature
{
  spectrafunc_ = spectrafunc;
}

void PlanckSpectra::SetSpectraVar(SpectraVar spectravar)
// spectravar = NU     : variable is frequency (Hz)
//            = LAMBDA : variable is wavelength (m)
{
  spectravar_ = spectravar;
}

void PlanckSpectra::SetSpectraPower(SpectraPower spectrapower)
// spectrapower = POWER : in W/...
//              = PHOTON : in photon/s/...
{
  spectrapower_ = spectrapower;
}

void PlanckSpectra::SetSpectraUnit(SpectraUnit spectraunit)
// spectraunit = ANGSFLUX  : radiance   W/m^2/Sr/Hz  or ...
//             = SFLUX     : emittance  W/m^2/Hz     or ... 
//             = DENSENERG : density    J/m^3/Hz     or ...
//             = EXPON     : just the function exp() or equivalent for Rayleigh
{
  spectraunit_ = spectraunit;
}


//************************************************************************
double PlanckSpectra::F2L(double nu_or_lambda)
// Convert frequency (GHz) to/from Wavelength (m)
// Lambda(m) * Nu(GHz) = C(km/s) * 10^-6
{
  return CinM_ / nu_or_lambda;
}

//************************************************************************
double PlanckSpectra::PlanckExp(double fl)
{
  double x = (spectravar_==LAMBDA) ? HCsK_/(fl*T_) : HsK_*fl/T_;

  if(spectraapprox_ == RAYLEIGH) return 1./x;  // Rayleigh

  if(x>MAXARGEXP) return 0.;

  if(spectraapprox_ == WIEN) return exp(-x);  // Wien

  if(x<MINARGEXP) return 1./(x*(1.+x/2.));
  return 1./(exp(x)-1.);
}

double PlanckSpectra::DPlanckExp_DT(double fl)
{
  if(spectraapprox_ == RAYLEIGH)
    return (spectravar_==LAMBDA) ? fl/HCsK_: 1. / (HsK_*fl);  // Rayleigh

  double x = (spectravar_==LAMBDA) ? HCsK_/(fl*T_) : HsK_*fl/T_;

  if(x>MAXARGEXP) return 0.;

  if(spectraapprox_ == WIEN) return x/T_ * exp(-x);  // Wien

  double v;
  if(x<MINARGEXP) { // 1/(exp(x)-1)  -derive-> -exp(x)/(exp(x)-1)^2
    v = x*(1.+x/2.);
    return x/T_ * (1.+v)/(v*v);
  }
  v = exp(x);
  return x/T_ * v/((v-1.)*(v-1.));  
}

//------------------------------------------------------------
double PlanckSpectra::operator() (double x)
{
  // --- type de spectre et d'approximation selon choix variable
  double fex = (spectrafunc_==DERIV) ? DPlanckExp_DT(x): PlanckExp(x);

  // --- le coefficient multiplicateur devant l'expo
  if(spectraunit_ == EXPON) return fex;  // expo seule

  // radiance en ph/...               longeur d'onde         frequency
  double coeff = (spectravar_==LAMBDA) ? 2.*CinM_/(x*x*x*x): 2.*x*x/(CinM_*CinM_);

  if(spectraunit_ == SFLUX)      coeff *= M_PI;  // emittance en ph/...
  else if(spectraunit_ == DENSENERG) coeff *= 4.*M_PI/CinM_;  // densite en ph/...

  // unite en W/...                longeur d'onde         frequency
  if(spectrapower_ == POWER)
    coeff *= (spectravar_==LAMBDA) ? h_Planck_Cst*CinM_/x: h_Planck_Cst*x;

  return coeff * fex;
}

//------------------------------------------------------------
double PlanckSpectra::PlanckEnergie(void)
//  densite d'energie: energie par m^3  en "J / m^3"
// cste de Stefan-Boltzmann en W/m^2/K^4
{
  return 4.*StefBoltz_Cst*T_*T_*T_*T_/(SpeedOfLight_Cst*1000.);
}

//------------------------------------------------------------
double PlanckSpectra::PlanckPhoton(void)
//  densite de photons: nombre de photons par m^3  en "ph / m^3"
{
  double x = k_Boltzman_Cst*T_/(h_Planck_Cst*SpeedOfLight_Cst*1000.);
  return 16.*M_PI*Zeta_3_Cst* x*x*x;
}

double PlanckSpectra::WienLaw(void)
// Retourne la frequence (ou longeur d'onde) du maximum
// du spectre de radiance (ou emittance ou densite).
// La valeur du maximum depend de l'unite du spectre (W/... ou Ph/s/...)
/*
--- Les fonctions:
En frequence: n = 3  en W/...  ,  n = 2  en Ph/...
En longueur d'onde: n = 5  en W/...  ,  n = 4  en Ph/...
f[x_] := x^n*Exp[-x]; Solve[f'[x]==0,x]   // Wien
f[x_] := x^n/(Exp[x]-1); NSolve[f'[x]==0,x,17]  // Planck
La derivee par rapport a T des fonctions en frequence:
df/dT[x_] := x^(n+1)*Exp[-x]; Solve[f'[x]==0,x]  // Wien
df/dT[x_] := x^(n+1)*Exp[x]/(Exp[x]-1)^2; NSolve[f'[x]==0,x,17]  // Planck
*/
{
  double x=-1.;
  if(spectravar_ == NU) { // Frequence
    if(spectrafunc_ == VALUE) {  // Le spectre
      if(spectrapower_ == POWER) //  W/...
        {if(spectraapprox_==PLANCK) x=2.82143938074; else if(spectraapprox_==WIEN) x=3.;}
      else if(spectrapower_ == PHOTON) //  Ph/s/...
        {if(spectraapprox_==PLANCK) x=1.59362426; else if(spectraapprox_==WIEN) x=2.;}
    } else if(spectrafunc_ == DERIV) {  // La derivee du spectre
      if(spectrapower_ == POWER) //  W/... 
        {if(spectraapprox_==PLANCK) x=3.830016095; else if(spectraapprox_==WIEN) x=4.;}
      else if(spectrapower_ == PHOTON) //  Ph/s/...
        {if(spectraapprox_==PLANCK) x=2.575678888; else if(spectraapprox_==WIEN) x=3.;}
    }
  } else { // Longeur d'onde
    if(spectrafunc_ == VALUE) {  // Le spectre
      if(spectrapower_ == POWER) //  W/...
        {if(spectraapprox_==PLANCK) x=4.96511420; else if(spectraapprox_==WIEN) x=5.;}
      else if(spectrapower_ == PHOTON) //  Ph/s/...
        {if(spectraapprox_==PLANCK) x=3.9206904099; else if(spectraapprox_==WIEN) x=4.;}
    } else if(spectrafunc_ == DERIV) {  // La derivee du spectre
      if(spectrapower_ == POWER) //  W/...
        {if(spectraapprox_==PLANCK) x=5.969409190; else if(spectraapprox_==WIEN) x=6.;}
     else if(spectrapower_ == PHOTON) //  Ph/s/...
        {if(spectraapprox_==PLANCK) x=4.928119345; else if(spectraapprox_==WIEN) x=5.;}
    }
  }

  if(x<0.) return x;  // Rayleigh n'a pas de maximum
  if(spectravar_ == NU) return x / HsK_ *T_; // Frequence
  return HCsK_ / x /T_; // Longeur d'onde
}

double PlanckSpectra::FindMaximum(double eps)
// Retourne le maximum des spectres (a la precicion eps)
{
  if(spectraapprox_ == RAYLEIGH) {  // Pas de maximum pour les spectres de Raleigh
    cout<<"PlanckSpectra::FindMaximum : no maximum for Raleigh spectra"<<endl;
    throw ParmError("PlanckSpectra::FindMaximum : no maximum for Raleigh spectra");
  }
  if(eps<=0.) eps = 0.001;

  double xmax = WienLaw(), wmax=(*this)(xmax);
  double xlim[2]={xmax,xmax};

  // Encadrement du maximum par xlim[0] et xlim[1]
  while( (*this)(xlim[1])>=0.01*wmax ) xlim[1] *= 2.;
  while( (*this)(xlim[0])>=0.01*wmax ) xlim[0] /= 2.;
  //cout<<"FindMaximum: wmax="<<wmax<<" pour x=["<<xlim[0]<<","<<xlim[1]<<"]"<<endl;

  while( (xlim[1]-xlim[0])/xmax > eps ) {
   wmax=-1.;
   double dx = (xlim[1]-xlim[0])/10.;
   for(double x=xlim[0]; x<=xlim[1]; x+=dx) {
     double f=(*this)(x);
     if(f>wmax) {wmax=f; xmax=x;}
   }
   //cout<<"...wmax="<<wmax<<" pour xmax="<<xmax<<" dx="<<dx<<endl;
   xlim[0] = xmax-dx*1.1; xlim[1] = xmax+dx*1.1;
 }


  //cout<<"wmax="<<wmax<<" pour x="<<xmax<<"  dx="<<xlim[1]-xlim[0]<<endl;
  return xmax;
}

}  // Fin namespace SOPHYA