#ifndef GRAYBODY_H
#define GRAYBODY_H

#include <cmath>

#include "Dust.h"
#include "manticore.h"

#ifndef CMB
#define CMB 0
#endif

namespace manticore {

// Forward references.
class Detector;


/**
  \brief %Graybody emission model.

  A graybody is a blackbody emitter modified by an optical depth factor:
\f[
  I^{\rm gray}_\nu = (1 - e^{-\tau})B_\nu(T),
\f]
  where \f$\tau = \tau(\nu; \Sigma)\f$ is the optical depth, \f$\Sigma\f$
  is the mass column density and \f$B_\nu(T)\f$ is the Planck function:
\f[
  B_\nu(T) = {2h\nu^3\over c^2}{1\over e^{h\nu/kT}-1}.
\f]

  The functional form for \f$\tau\f$ depends on the details of the dust model.
  This class accepts an arbitrary dust opacity model.

  Both single- and dual-temperature emission models are supported.
  Single-temperature models contain a simple uniform slab of material.
  Dual-temperature models consist of a slab of "cold" (or "core") material
  sandwiched between two identical layers of "hot" (or "halo") material.
*/
class Graybody {
public:
  /// Blackbody convenience value.
  static constexpr double B0 = 2.0*h_Planck/(c_light*c_light);

  /// CMB temperature (K).
  static constexpr double Tcmb = 2.73;

  /// Single dust model constructor.
  /// \param[in] dust Dust model.
  /// \param[in] err  Relative error target for integrations.
  Graybody(Dust &dust, double err = 1e-3) :
    dust_{dust, dust}, dualDust_{false}, err_{err},
    hacc_{gsl_interp_accel_alloc()}, cacc_{gsl_interp_accel_alloc()} { }

  /// Dual dust model constructor.
  /// \param[in] dusth Hot/halo dust model.
  /// \param[in] dustc Cold/core dust model.
  /// \param[in] err   Relative error target for integrations.
  Graybody(Dust &dusth, Dust &dustc, double err = 1e-3) :
    dust_{dusth, dustc}, dualDust_{dusth.name() != dustc.name()}, err_{err},
    hacc_{gsl_interp_accel_alloc()}, cacc_{gsl_interp_accel_alloc()} { }

  /// Copy constructor.
  Graybody(const Graybody &g) :
    dust_{g.dust_}, dualDust_{g.dualDust_}, err_{g.err_},
    hacc_{gsl_interp_accel_alloc()}, cacc_{gsl_interp_accel_alloc()} { }

  /// Destructor.
  ~Graybody() {
    gsl_interp_accel_free(hacc_);
    gsl_interp_accel_free(cacc_);
  }

  /// Blackbody specific intensity (erg/cm^2/s/Hz/sr).
  /// \param[in] nu  Frequency (Hz).
  /// \param[in] T   Temperature (K).
  double Bnu(double nu, double T) const
  { return B0*nu*nu*nu/expm1(h_Planck*nu/(k_Boltzman*T)); }

  /// Blackbody specific intensity T-derivative (erg/cm^2/s/Hz/sr/K).
  /// \param[in] nu  Frequency (Hz).
  /// \param[in] T   Temperature (K).
  double dBnu_dT(double nu, double T) const {
    double x = h_Planck*nu/(k_Boltzman*T), em1 = expm1(x);
    return B0*nu*nu*nu*(em1+1.0)*x/(em1*em1*T);
  }

  /// %Single-temperature graybody specific intensity (erg/cm^2/s/Hz/sr).
  /// \param[in] nu     Frequency (Hz).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  double Inu(double nu, double T, double Sigma) const {
    double em1 = expm1(-dust_.first.kappa(nu, hacc_)*Sigma);
    return -em1*(Bnu(nu, T)-(CMB ? Bnu(nu, Tcmb) : 0.0));
  }

  /// %Single-temperature graybody specific intensity T-derivative
  ///  (erg/cm^2/s/Hz/sr/K).
  /// \param[in] nu     Frequency (Hz).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  double dInu_dT(double nu, double T, double Sigma) const
  { return -expm1(-dust_.first.kappa(nu, hacc_)*Sigma)*dBnu_dT(nu, T); }

  /// %Single-temperature graybody specific intensity Sigma-derivative
  ///  (erg/s/Hz/sr/g).
  /// \param[in] nu     Frequency (Hz).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  double dInu_dS(double nu, double T, double Sigma) const {
    double kappa = dust_.first.kappa(nu, hacc_);
    return kappa*exp(-kappa*Sigma)*(Bnu(nu, T)-(CMB ? Bnu(nu, Tcmb) : 0.0));
  }

  /// %Two-temperature graybody specific intensity (erg/cm^2/s/Hz/sr).
  /// \param[in] nu      Frequency (Hz).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  double Inu(double nu, double Tc, double Sigmac,
                        double Th, double Sigmah) const
  {
    double kappah = dust_.first.kappa(nu, hacc_),
           kappac = (dualDust_ ? dust_.second.kappa(nu, cacc_) : kappah),
             ehm1 = expm1(-kappah*Sigmah), eh = 1.0 + ehm1,
             ecm1 = expm1(-kappac*Sigmac), ec = 1.0 + ecm1;
    return -ehm1*(1.0+eh*ec)*Bnu(nu, Th)-ecm1*eh*Bnu(nu, Tc);
  }

  /// %Two-temperature graybody specific intensity Tc-derivative
  ///  (erg/cm^2/s/Hz/sr/K).
  /// \param[in] nu      Frequency (Hz).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  double dInu_dTc(double nu, double Tc, double Sigmac,
                             double Th, double Sigmah) const {
    (void )Th; // UNUSED
    double kappah = dust_.first.kappa(nu, hacc_),
           kappac = (dualDust_ ? dust_.second.kappa(nu, cacc_) : kappah);
    return -expm1(-kappac*Sigmac)*exp(-kappah*Sigmah)*dBnu_dT(nu, Tc);
  }

  /// %Two-temperature graybody specific intensity Sigmac-derivative
  ///  (erg/s/Hz/sr/g).
  /// \param[in] nu      Frequency (Hz).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  double dInu_dSc(double nu, double Tc, double Sigmac,
                             double Th, double Sigmah) const {
    double kappah = dust_.first.kappa(nu, hacc_),
           kappac = (dualDust_ ? dust_.second.kappa(nu, cacc_) : kappah);
    return kappac*exp(-(kappah*Sigmah+kappac*Sigmac)) *
           (expm1(-kappah*Sigmah)*Bnu(nu, Th) + Bnu(nu, Tc));
  }

  /// %Two-temperature graybody specific intensity Th-derivative
  ///  (erg/cm^2/s/Hz/sr/K).
  /// \param[in] nu      Frequency (Hz).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  double dInu_dTh(double nu, double Tc, double Sigmac,
                             double Th, double Sigmah) const {
    (void )Tc; // UNUSED
    double kappah = dust_.first.kappa(nu, hacc_),
           kappac = (dualDust_ ? dust_.second.kappa(nu, cacc_) : kappah),
           ehc = exp(-(kappah*Sigmah+kappac*Sigmac));
    return -expm1(-kappah*Sigmah)*(1.0+ehc)*dBnu_dT(nu, Th);
  }

  /// %Two-temperature graybody specific intensity Sigmah-derivative
  ///  (erg/s/Hz/sr/g).
  /// \param[in] nu      Frequency (Hz).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  double dInu_dSh(double nu, double Tc, double Sigmac,
                             double Th, double Sigmah) const {
    double kappah = dust_.first.kappa(nu, hacc_),
           kappac = (dualDust_ ? dust_.second.kappa(nu, cacc_) : kappah),
             ehm1 = expm1(-kappah*Sigmah), eh = 1.0 + ehm1,
             ecm1 = expm1(-kappac*Sigmac), ec = 1.0 + ecm1;
    return kappah*eh*((2.0*eh*ec-ecm1)*Bnu(nu, Th)+ecm1*Bnu(nu, Tc));
  }

  /// %Single-temperature graybody integrated flux (erg/cm^2/s).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  /// \param[in] detect Detector instance.
  /// \param[in] sr     Solid angle (sr, optional).
  double F(double T, double Sigma,
           const Detector *detect, double sr = 1.0) const
  { return Itot(&Graybody::Inu, T, Sigma, detect, sr); }

  /// %Single-temperature graybody integrated flux T-derivative (erg/cm^2/s/K).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  /// \param[in] detect Detector instance.
  /// \param[in] sr     Solid angle (sr, optional).
  double dF_dT(double T, double Sigma,
               const Detector *detect, double sr = 1.0) const
  { return Itot(&Graybody::dInu_dT, T, Sigma, detect, sr); }

  /// %Single-temperature graybody integrated flux Sigma-derivative (erg/s/g).
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  /// \param[in] detect Detector instance.
  /// \param[in] sr     Solid angle (sr, optional).
  double dF_dS(double T, double Sigma,
               const Detector *detect, double sr = 1.0) const
  { return Itot(&Graybody::dInu_dS, T, Sigma, detect, sr); }

  /// %Two-temperature graybody integrated flux (erg/cm^2/s).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double F(double Tc, double Sigmac, double Th, double Sigmah,
           const Detector *detect, double sr = 1.0) const
  { return Itot2(&Graybody::Inu, Tc, Sigmac, Th, Sigmah, detect, sr); }

  /// %Two-temperature graybody integrated flux Tc-derivative (erg/cm^2/s/K).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double dF_dTc(double Tc, double Sigmac, double Th, double Sigmah,
                const Detector *detect, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dTc, Tc, Sigmac, Th, Sigmah, detect, sr); }

  /// %Two-temperature graybody integrated flux Sigmac-derivative (erg/s/g).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double dF_dSc(double Tc, double Sigmac, double Th, double Sigmah,
                const Detector *detect, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dSc, Tc, Sigmac, Th, Sigmah, detect, sr); }

  /// %Two-temperature graybody integrated flux Th-derivative (erg/cm^2/s/K).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double dF_dTh(double Tc, double Sigmac, double Th, double Sigmah,
                const Detector *detect, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dTh, Tc, Sigmac, Th, Sigmah, detect, sr); }

  /// %Two-temperature graybody integrated flux Sigmah-derivative (erg/s/g).
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double dF_dSh(double Tc, double Sigmac, double Th, double Sigmah,
                const Detector *detect, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dSh, Tc, Sigmac, Th, Sigmah, detect, sr); }

protected:
  /// %Single-temperature graybody integrated flux variant (erg/cm^2/s/??).
  /// \param[in] I      Specific flux integration function.
  /// \param[in] T      Temperature (K).
  /// \param[in] Sigma  Mass density (g/cm^2).
  /// \param[in] detect Detector instance.
  /// \param[in] sr     Solid angle (sr, optional).
  double Itot(double (Graybody::*I)(double, double, double) const,
              double T, double Sigma,
              const Detector *detect, double sr = 1.0) const;

  /// %Two-temperature graybody integrated flux variant (erg/cm^2/s/??).
  /// \param[in] I       Specific flux integration function.
  /// \param[in] Tc      Cold/core temperature (K).
  /// \param[in] Sigmac  Cold/core mass density (g/cm^2).
  /// \param[in] Th      Hot/halo  temperature (K).
  /// \param[in] Sigmah  Hot/halo  mass density (g/cm^2).
  /// \param[in] detect  Detector instance.
  /// \param[in] sr      Solid angle (sr, optional).
  double Itot2(double (Graybody::*I)(double,double,double,double,double) const,
               double Tc, double Sigmac, double Th, double Sigmah,
               const Detector *detect, double sr = 1.0) const;

  /// Dust models.
  std::pair<Dust&, Dust&> dust_;

  /// Whether the two dust models differ;
  bool dualDust_;

  /// Integration relative error.
  double err_;

  /// Interpolation index accelerators.
  gsl_interp_accel
    *hacc_ = nullptr,  //< Dust::kappa() lookup accelerator (halo).
    *cacc_ = nullptr;  //< Dust::kappa() lookup accelerator (core).
};

}  // namespace manticore

#endif  // GRAYBODY_H