manticore.cc

Go to the documentation of this file.

#include "Detector.h"
#include "Graybody.h"
#include "PACS.h"
#include "version.h"

namespace mu = mutils;

namespace manticore {

void ccfCheck()
{
  auto f = mu::Log::stream;
  bool extend = false;
  const double nuHz_1micron = 1e4*c_light;
  Dust dust(1.0, 1e10, 0.0);
  Graybody gray(dust, 1e-4);  // Blackbody.
  PACS pacs70(70, extend), pacs100(100, extend), pacs160(160, extend);
  pacs70.init();
  pacs100.init();
  pacs160.init();

  fputs("\n**** Effective bandwidth (Hz) ****\n", f);
  fputs("   70u     100u      160u\n", f);
  fprintf(f, "%8.2e  %8.2e  %8.2e\n",
          pacs70.freqWidth(), pacs100.freqWidth(), pacs160.freqWidth());

  fputs("\n**** Blackbody color corrections ****\n", f);
  fputs("Temp(K)    70u    100u    160u\n", f);
  for (auto T: {1e4, 1e3, 1e2, 5e1, 1e1, 5e0}) {
    fprintf(f, "%5gK %7.3f %7.3f %7.3f\n", T,
            pacs70.ccf(gray, T), pacs100.ccf(gray, T), pacs160.ccf(gray, T));
  }

  fputs("\n**** Graybody color corrections ****\n", f);
  fputs("(beta, T)     70u  100u  160u\n", f);
  for (double beta=1.0; beta <= 2.0; beta++) {
    dust.setModel(nuHz_1micron/1.0, 1.0, beta);
    for (auto T: {1e2, 1e1}) {
      fprintf(f, "(%g, %3gK) = %.3f %.3f %.3f\n", beta, T,
              pacs70.ccf(gray, T), pacs100.ccf(gray, T), pacs160.ccf(gray, T));
    }
  }

  double T = 15.0, Sigma = 1.0, sr = 3.62e-9;
  dust.setModel(nuHz_1micron/157.0, 1.0, 1.5);
  fprintf(f, "\nExample flux = %.2e erg/cm^2/s, Fnu = %.2f Jy, 1/CCF = %.3f\n",
          gray.F(T, Sigma, &pacs100, sr),
          gray.Inu(pacs100.freqBand(), T, Sigma)*sr/1e-23,
          1.0/pacs100.ccf(gray, T, Sigma));

  fputs("\n**** Derivative test ****\n", f);
  fprintf(f, "dF/dT: %.4e calc, %.4e true\n",
          (gray.F(T+1e-3, Sigma, &pacs100, sr)-
           gray.F(T-1e-3, Sigma, &pacs100, sr))/2e-3,
           gray.dF_dT(T, Sigma, &pacs100, sr));
  fprintf(f, "dF/dS: %.4e calc, %.4e true\n",
          (gray.F(T, Sigma+1e-3, &pacs100, sr)-
           gray.F(T, Sigma-1e-3, &pacs100, sr))/2e-3,
           gray.dF_dS(T, Sigma, &pacs100, sr));
}

gsl_spline *initSpline(const tableType &table, bool ln, unsigned stride)
{
  std::vector<double> x, y;
  for (unsigned i=0; i < table.size(); i += stride) {
    auto &xy = table[i];
    x.push_back(xy.first), y.push_back(ln ? log(xy.second) : xy.second);
  }
  auto spline = gsl_spline_alloc(gsl_interp_akima, x.size());
  gsl_spline_init(spline, &x[0], &y[0], x.size());
  return spline;
}

void printVersion()
{
  const char *what;
  std::string when;
  if (*manticore::reldate) { what = "released", when = manticore::reldate; }
  else { what = "built", when = std::string(__DATE__) + " " + __TIME__; }
  fprintf(mu::Log::stream, "Manticore version %s, %s %s.\n\n",
          manticore::version, what, when.c_str());
}

int getBand(const CCfits::HDU &hdu)
{
  double band;
  return (int )hdu.keyWord("BAND").value(band);
}

}  // namespace manticore


int main(int argc, char *argv[])
try {
  // Command line interface.
  mu::CommandLine cli;
  std::vector<mu::CommandLine::Option> opts{
    {'1', "one-param", "", "", "", "",
     "Fit one parameter dust emission model (default: two-param).\n"
     "(Fixed-temperature cloud model.)"},

    {'2', "two-param", "", "", "", "",
     "Fit two parameter dust emission model (default).\n"
     "(Single-temperature cloud model.)"},

    {'3', "three-param", "", "", "", "",
     "Fit three parameter dust emission model (default: two-param).\n"
     "(Two-temperature cloud model with prescribed halo temperature.)"},

    {'4', "four-param", "", "", "", "",
     "Fit four parameter dust emission model (default: two-param).\n"
     "(Two-temperature cloud model, all parameters free.)"},

    {'A', "accuracy", "RELERR", "", "", "",
     "Maximum relative error for least-squares fit (default: 1e-3)."},

    {'a', "abs-diff", "", "", "", "",
     "Output absolute differences in delta maps (default: sigma-relative)."},

    {'B', "no-band-output", "", "", "", "",
     "Do not append input band data to output FITS file."},

    {'C', "cutout", "COLSxROWS@COL0,ROW0", "", "", "",
     "Process only a COLSxROWS cutout at lower-left corner (COL0,ROW0)."},

    {'c', "chi2-max", "CHI2", "", "", "",
     "Maximum chi^2/DOF for meaningful 2-param result (default: 1.0)."},

    {'D', "dust-model", "(MODEL | MODELhalo,MODELcore)", "", "", "",
     "Dust opacity model(s) (default: DL3).\n"
     "Use pow:nu0_GHz:kappa0:beta for power law kappa0*(nu/nu0)^beta."},

    {'e', "error-map", "FILE", "", "", "",
     "Input FITS error map (BAND and TELESCOP keys must match parent image)."},

    {'E', "error-pct", "PCT", "", "", "",
     "Percent flux error when error map unavailable (default: 3.0)."},

    {'F', "force", "", "", "", "",
     "Force overwrite of existing files."},

    {'G', "grid-samp", "", "", "", "",
     "Sample data on a fixed grid to estimate parameter variance.\n"
     "(Total number of samples is 2^numBands.)"},

    {'g', "gas-to-dust", "(RATIO | RATIOhalo,RATIOcore)", "", "", "",
     "Gas-to-dust ratio(s) (default: 100.0)."},

    {'h', "help", "", "", "", "",
     "Display command line help summary and exit."},

    {'l', "log-file", "FILE", "", "", "",
     "Redirect log output to FILE (default is stderr)."},

    {'M', "model-map", "", "", "", "",
     "Replace input data with model observation grid."},

    {'m', "max-iterations", "NMAX", "", "", "",
     "Maximum least-squares iterations per pixel (default: 50)."},

    {'N', "nh2-min", "NH2MIN", "", "", "",
     "Minimum allowed (halo) N(H2) column density in cm^-2."},

    {'n', "num-bands", "NMIN", "", "", "",
     "Minimum number of bands required for least-squares fit (default: 3)."},

    {'o', "output-file", "FILE", "", "", "",
     "Output result to FITS file FILE (default: mtcore.fits)."},

    {'P', "pacs-check", "", "", "", "",
     "Generate comparative PACS color correction tables and exit."},

    {'p', "point-source", "", "", "", "",
     "Use point source detector response (default is extended sources)."},

    {'q', "quiet", "", "", "", "",
     "Quiet mode; decrement logging level (opposite of -v).\n"
     "(Specify multiple times for added effect.)"},

    {'R', "rand-samp", "N", "", "", "",
     "Randomly sample data N times to estimate parameter variance.\n"
     "(Without -G, -R or -S, the covariance matrix is used instead.)"},

    {'S', "sys-samp", "", "", "", "",
     "Sample systematic data errors to estimate parameter variance.\n"
     "(May be combined with -G or -R.)"},

    {'s', "single-temp", "FILE", "", "", "",
     "FITS output file from the single-temperature (2-param) fit."},

    {'T', "temp0", "T0", "", "", "",
     "Initial guess for (halo) temperature in Kelvin (default: 15.0)."},

    {'t', "threads", "NTHREADS", "", "", "",
     "Use NTHREADS threads for processing (default: 1)."},

    {'v', "verbose", "", "", "", "",
     "Verbose mode; increment logging level (opposite of -q).\n"
     "(Specify multiple times for added effect.)"},

    {'V', "version", "", "", "", "",
     "Display program version and exit."},
  };

  // Parse command line options.
  cli.init(argc, argv, &opts);
  manticore::dust::printModels(opts, 'D');

  // --quiet / --verbose
  int verb = cli.getOption('v') - cli.getOption('q');
  mu::Log::level += verb;
  CCfits::FITS::setVerboseMode(verb > 0);

  // --help
  if (cli.getOption('h')) {
    bool quiet = (mu::Log::level <= mu::Log::WARN);
    if (!quiet) { manticore::printVersion(); }
    cli.usage(mu::Log::stream, opts, "BAND1.fits BAND2.fits BAND3.fits [...]",
              quiet);
    if (quiet) {
      fputs("\nFor more details, use --verbose --help.\n", mu::Log::stream);
    }
    exit(EXIT_SUCCESS);
  };

  // --version
  if (cli.getOption('V')) {
    manticore::printVersion();
    exit(EXIT_SUCCESS);
  }

  // --log-file
  mu::Log::setStream(cli.getString('l', "!"));

  // --check
  if (cli.getOption('P')) {
    manticore::ccfCheck();
    exit(EXIT_SUCCESS);
  }

  manticore::process(cli);

  return EXIT_SUCCESS;
} catch (std::exception &ex) {
  mu::printException(ex);
  return EXIT_FAILURE;
}