Graybody.cc

Go to the documentation of this file.

#include "mutils/nmath.h"

#include "Detector.h"
#include "Graybody.h"

namespace mu = mutils;

namespace manticore {

namespace graybody {

extern "C" {

struct FluxData {
  double Tc,      
         Sigmac,  
         Th,      
         Sigmah;  

  const Detector *d;  

  const Graybody *g;  

  double (Graybody::*I)(double, double, double) const;

  double (Graybody::*I2)(double, double, double, double, double) const;
};

double flux(double nu, void *fd) {
  auto &p = *static_cast<const FluxData*>(fd);
  return (p.I2 ? (p.g->*(p.I2))(nu, p.Tc, p.Sigmac, p.Th, p.Sigmah) :
                 (p.g->*(p.I ))(nu, p.Tc, p.Sigmac)) *
         (p.d ? p.d->response(nu) : 1.0);
}

}  // extern "C"

}  // namescape graybody


double Graybody::Itot(double (Graybody::*I)(double, double, double) const,
                      double T, double Sigma, const Detector *detect,
                      double sr) const
{
  std::pair<double, double> freqRange =
    (detect ? detect->freqRange() :
              std::pair<double,double>{0.0,12.0*k_Boltzman*T/h_Planck});

  double result, abserr;
  graybody::FluxData fd = {T, Sigma, 0.0, 0.0, detect, this, I, nullptr};
  gsl_function f = {graybody::flux, &fd};
  gsl_integration_qag(&f, freqRange.first, freqRange.second,
                      0.0, err_, subLimit, intKey,
                      integ_, &result, &abserr);
  return result*sr;
}


double Graybody::Itot2(
         double (Graybody::*I)(double, double, double, double, double) const,
         double Tc, double Sigmac, double Th, double Sigmah,
         const Detector *detect, double sr) const
{
  std::pair<double, double> freqRange =
    (detect ? detect->freqRange() :
              std::pair<double,double>{0.0, 12.0*k_Boltzman*std::max(Tc, Th)/
                                                 h_Planck});

  graybody::FluxData fd = {Tc, Sigmac, Th, Sigmah, detect, this, nullptr, I};
  gsl_function f = {graybody::flux, &fd};
  /*
  double result, abserr;
  gsl_integration_qag(&f, freqRange.first, freqRange.second,
                      0.0, err_, subLimit, intKey,
                      integ_, &result, &abserr);
  return sr*result;
  */
//  MU_WARN_IF(Tc < 1.0, "Graybody::Itot2", "Low Tc = %.4f", Tc);
  return sr * mu::kgauss_r(f.function, f.params,
                           freqRange.first, freqRange.second, 0.1*err_);
}

}  // namespace manticore