Graybody.h

Go to the documentation of this file.

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


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

  static constexpr double Tcmb = 2.73;

  Graybody(Dust &dust, double err = 1e-3) :
    dust_{dust, dust}, dualDust_{false}, err_{err},
    hacc_{gsl_interp_accel_alloc()}, cacc_{gsl_interp_accel_alloc()} { }

  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()} { }

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

  ~Graybody() {
    gsl_interp_accel_free(hacc_);
    gsl_interp_accel_free(cacc_);
  }

  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 {
    double em1 = expm1(-dust_.first.kappa(nu, hacc_)*Sigma);
    return -em1*(Bnu(nu, T)-(CMB ? Bnu(nu, Tcmb) : 0.0));
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  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:
  double Itot(double (Graybody::*I)(double, double, double) const,
              double T, double Sigma,
              const Detector *detect, double sr = 1.0) const;

  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;

  std::pair<Dust&, Dust&> dust_;

  bool dualDust_;

  double err_;

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

}  // namespace manticore

#endif  // GRAYBODY_H