#ifndef GRAYBODY_H
#define GRAYBODY_H

#include <gsl/gsl_integration.h>

#include <cmath>

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

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.
*/
class Graybody {
public:
  /// Blackbody convenience value.
  static constexpr double B0 = 2.0*h_Planck/(c_light*c_light);

  /// Default constructor.
  /// \param[in] dust Dust model.
  /// \param[in] err  Relative error target for integrations.
  Graybody(Dust &dust, double err = 1e-3) :
    dust_(dust), err_(err), integ_(gsl_integration_workspace_alloc(subLimit))
    { }

  /// Copy constructor.
  Graybody(const Graybody &g) : dust_(g.dust_), err_(g.err_),
    integ_(gsl_integration_workspace_alloc(subLimit)) { }

  /// Dust model reference.
  Dust &dust() noexcept { return dust_; }

  /// Dust model reference (const).
  const Dust &dust() const noexcept { return dust_; }

  /// 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);
  }

  /// %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
  { return -expm1(-dust_.kappa(nu)*Sigma)*Bnu(nu, T); }

  /// %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_.kappa(nu)*Sigma)*dBnu_dT(nu, T); }

  /// %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_.kappa(nu);
    return exp(-kappa*Sigma)*kappa*Bnu(nu, T);
  }

  /// %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 (optional).
  /// \param[in] sr     Solid angle (sr, optional).
  double F(double T, double Sigma,
           const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot(&Graybody::Inu, T, Sigma, detect, sr); }

  /// %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 (optional).
  /// \param[in] sr     Solid angle (sr, optional).
  double dF_dT(double T, double Sigma,
               const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot(&Graybody::dInu_dT, T, Sigma, detect, sr); }

  /// %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 (optional).
  /// \param[in] sr     Solid angle (sr, optional).
  double dF_dS(double T, double Sigma,
               const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot(&Graybody::dInu_dS, T, Sigma, detect, sr); }

protected:
  /// Integral subdivision limit.
  static constexpr size_t subLimit = 512;

  /// Integration rule key.
  static constexpr int intKey = GSL_INTEG_GAUSS41;

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

  /// Dust model.
  Dust &dust_;

  /// Integration relative error.
  double err_;

  /// Integration workspace.
  gsl_integration_workspace *integ_;
};

}  // namespace manticore

#endif  // GRAYBODY_H