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

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;

  Dust &dust_;

  double err_;

  gsl_integration_workspace *integ_;
};

}  // namespace manticore

#endif  // GRAYBODY_H