Graybody.h

Go to the documentation of this file.

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


class Graybody {
public:
  static constexpr double B0 = 2.0*h_Planck/(c_light*c_light);

  Graybody(Dust &dust, double err = 1e-3) :
    dust_(dust), err_(err), integ_(gsl_integration_workspace_alloc(subLimit))
    { }

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

  Dust &dust() noexcept { return dust_; }

  const Dust &dust() const noexcept { return dust_; }

  double Bnu(double nu, double T) const
  { return B0*nu*nu*nu/expm1(h_Planck*nu/(k_Boltzman*T)); }

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

  double Inu(double nu, double T, double Sigma) const
  { return -expm1(-dust_.kappa(nu)*Sigma)*Bnu(nu, T); }

  double dInu_dT(double nu, double T, double Sigma) const
  { return -expm1(-dust_.kappa(nu)*Sigma)*dBnu_dT(nu, T); }

  double dInu_dS(double nu, double T, double Sigma) const {
    double kappa = dust_.kappa(nu);
    return exp(-kappa*Sigma)*kappa*Bnu(nu, T);
  }

  double Inu(double nu, double Tc, double Sigmac,
                        double Th, double Sigmah) const
  {
    double kappa = dust_.kappa(nu),
            ecm1 = expm1(-kappa*Sigmac), ec = 1.0 + ecm1,
            ehm1 = expm1(-kappa*Sigmah), eh = 1.0 + ehm1;
    return -ehm1*(1.0+eh*ec)*Bnu(nu, Th)-ecm1*eh*Bnu(nu, Tc);
  }

  double dInu_dTc(double nu, double Tc, double Sigmac,
                             double Th = 0.0, double Sigmah = 0.0) const {
    (void )Th; // UNUSED
    double kappa = dust_.kappa(nu);
    return -expm1(-kappa*Sigmac)*exp(-kappa*Sigmah)*dBnu_dT(nu, Tc);
  }

  double dInu_dSc(double nu, double Tc, double Sigmac,
                             double Th = 0.0, double Sigmah = 0.0) const {
    double kappa = dust_.kappa(nu);
    return kappa*exp(-kappa*(Sigmah+Sigmac)) *
           (expm1(-kappa*Sigmah)*Bnu(nu, Th) + Bnu(nu, Tc));
  }

  double dInu_dTh(double nu, double Tc, double Sigmac,
                             double Th = 0.0, double Sigmah = 0.0) const {
    (void )Tc; // UNUSED
    double kappa = dust_.kappa(nu), ehc = exp(-kappa*(Sigmah+Sigmac));
    return -expm1(-kappa*Sigmah)*(1.0+ehc)*dBnu_dT(nu, Th);
  }

  double dInu_dSh(double nu, double Tc, double Sigmac,
                             double Th = 0.0, double Sigmah = 0.0) const {
    double kappa = dust_.kappa(nu),
            ecm1 = expm1(-kappa*Sigmac), ec = 1.0 + ecm1,
            ehm1 = expm1(-kappa*Sigmah), eh = 1.0 + ehm1;
    return kappa*eh*((2.0*eh*ec-ecm1)*Bnu(nu, Th)+ecm1*Bnu(nu, Tc));
  }

  double F(double T, double Sigma,
           const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot(&Graybody::Inu, T, Sigma, detect, sr); }

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

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

  double F(double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
           const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot2(&Graybody::Inu, Tc, Sigmac, Th, Sigmah, detect, sr); }

  double dF_dTc(double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
                const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dTc, Tc, Sigmac, Th, Sigmah, detect, sr); }

  double dF_dSc(double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
                const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dSc, Tc, Sigmac, Th, Sigmah, detect, sr); }

  double dF_dTh(double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
                const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dTh, Tc, Sigmac, Th, Sigmah, detect, sr); }

  double dF_dSh(double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
                const Detector *detect = nullptr, double sr = 1.0) const
  { return Itot2(&Graybody::dInu_dSh, Tc, Sigmac, Th, Sigmah, detect, sr); }

protected:
  static constexpr size_t subLimit = 512;

  static constexpr int intKey = GSL_INTEG_GAUSS41;

  double Itot(double (Graybody::*I)(double, double, double) const,
              double T, double Sigma,
              const Detector *detect = nullptr, double sr = 1.0) const;

  double Itot2(double (Graybody::*I)(double,double,double,double,double) const,
               double Tc, double Sigmac, double Th = 0.0, double Sigmah = 0.0,
               const Detector *detect = nullptr, double sr = 1.0) const;

  Dust &dust_;

  double err_;

  gsl_integration_workspace *integ_;
};

}  // namespace manticore

#endif  // GRAYBODY_H