pacsResponse.cc

Go to the documentation of this file.

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