/**
  \file

  \brief Detector characteristics.

  \who \kpr
*/
#include "Detector.h"
#include "Graybody.h"

namespace mu = mutils;

namespace manticore {

/// Detector support.
namespace detector {

extern "C" {

/// Reference bandwidth integrand.
/// \param[in] nu Frequency (Hz).
/// \param[in] d  Detector instance.
double fWidth(double nu, void *d)
{ return static_cast<Detector*>(d)->response(nu)/nu; }

}

}  // namescape detector


void Detector::init()
{
  const char *const fn = "Detector::init";

  // Effective band limits.
  freqRange_ = {resp_.front().first, resp_.back().first};
  MU_DEBUG(fn, "%s response range is [%.2f, %.2f] microns", name(),
           1e4*c_light/freqRange_.second, 1e4*c_light/freqRange_.first);

  // Initialize integrator workspace.
  integ_ = gsl_integration_workspace_alloc(subLimit);

  // Initialize response interpolator.
  spline_ = initSpline(resp_);

  // Compute and set reference bandwidth.
  calcWidth();
}


/**
  \note
  For extended source detection, this response includes a correction for the
  frequency-dependent beam size (scaled to unity at freqBand()). The default
  implementation is based on linear interpolation of a pre-calculated response
  table.
*/
double Detector::response(double nu) const
{
  if (nu >= freqRange_.first && nu <= freqRange_.second) {
    return gsl_spline_eval(spline_, nu, nullptr);
  } else {
    // Truncate respone at specified band edges.
    return 0.0;
  }
}


/**
  Uses the response() function to compute the reference bandwidth (see detailed
  description).

  \post
  The #width_ member is set to the result.
*/
void Detector::calcWidth()
{
  const char *const fn = "Detector::calcWidth";
  double result, abserr;
  gsl_function f = {detector::fWidth, this};
  gsl_integration_qag(&f, freqRange_.first, freqRange_.second,
                      0.0, 1e-4, subLimit, intKey,
                      integ_, &result, &abserr);

  width_ = result*center_;
  MU_DEBUG(fn, "%s reference bandwidth: %.2f microns, %.3e Hz",
           name(), width_*1e4*c_light/(center_*center_), width_);
}


/**
  \note
  You must *divide* the reported detector flux density by the correction
  factor to obtain a fair estimate of the model flux density.
*/
double Detector::ccf(const Graybody &gray, double T, double Sigma) const
{ return gray.F(T, Sigma, this)/(width_*gray.Inu(center_, T, Sigma)); }

}  // namespace manticore