/**
  \file

  \brief Computes net response of a PACS band.

  The PACS bolometer response curve provided by Herschel support is
  provided over a completely different set of points than the supplied
  filter response tables. This program interpolates the bolometer response
  onto the filter grid and outputs the net response.

  \warning
  The input files must be 2-column {lambda, response} values only, no headers.
  The output is {nu, response} formatted for inclusion to Detector::tableType.

  \who \kpr
*/
#include <gsl/gsl_interp.h>
#include <mutils/util.h>
#include <manticore.h>

namespace mu = mutils;

int main(int argc, char *argv[])
{
  if (argc != 3) {
    fprintf(stderr, "Usage: %s PacsAbsorption.txt PacsFilter.txt\n", argv[0]);
    exit(EXIT_FAILURE);
  }

  // Open files.
  FILE *bolo = fopen(argv[1], "r"), *band = fopen(argv[2], "r");
  MU_EXCEPTION_ERRNO_IF(!bolo, argv[0], "Failed to open '%s'", argv[1]);
  MU_EXCEPTION_ERRNO_IF(!band, argv[0], "Failed to open '%s'", argv[2]);

  // Read tables.
  std::vector<double> lBolo, rBolo, lBand, rBand;
  double lamb, resp;
  while (fscanf(bolo, "%lf%lf", &lamb, &resp) == 2) {
    lBolo.push_back(lamb), rBolo.push_back(resp);
  }
  while (fscanf(band, "%lf%lf", &lamb, &resp) == 2) {
    lBand.push_back(lamb), rBand.push_back(resp);
  }
  fclose(bolo);
  fclose(band);

  // Interpolate and multiply.
  gsl_interp_accel *acc = gsl_interp_accel_alloc();
  gsl_interp *interp = gsl_interp_alloc(gsl_interp_cspline, lBolo.size());
  MU_EXCEPTION_IF(gsl_interp_init(interp, &lBolo[0], &rBolo[0], lBolo.size()),
                  argv[0], "Interpolator initialization failed");

  puts("  // {nu0, response0}");
  for (size_t i = lBand.size()-1; i != 0; i--) {
    double nu = 1e4*manticore::c_light / lBand[i];
    if (rBand[i] > 1e-4) {
      resp = rBand[i] *
             gsl_interp_eval(interp, &lBolo[0], &rBolo[0], lBand[i], acc);
      if (resp > 1e-4) { printf("  {%.5e, %.5e},\n", nu, resp); }
    }
  }

  exit(EXIT_SUCCESS);
}