/**
  \file

  \brief Dust mass and temperature fitting program.

  \who \kpr
*/
#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;
  Dust dust(1.0, 1e10, 0.0);
  Graybody gray(dust);  // Blackbody.
  PACS pacs70(70, extend), pacs100(100, extend), pacs160(160, extend);
  pacs70.init();
  pacs100.init();
  pacs160.init();

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

/**
  \brief Initializes cubic spline interpolator.

  \note
  The input \p table is not required for interpolation following
  initialization.

  \param[in] table  Interpolation table.
  \param[in] ln     Whether to take natural logarithm of the table ordinates.
  \return Initialized GSL spline instance.
*/
gsl_spline *initSpline(const tableType &table, bool ln)
{
  std::vector<double> x, y;
  for (auto &xy: table) {
    x.push_back(xy.first), y.push_back(ln ? log(xy.second) : xy.second);
  }
  auto spline = gsl_spline_alloc(gsl_interp_cspline, table.size());
  gsl_spline_init(spline, &x[0], &y[0], x.size());
  return spline;
}

/// Prints version information to the log.
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());
}

}  // namespace manticore


int main(int argc, char *argv[])
try {
  // Command line interface.
  mu::CommandLine cli;
  std::vector<mu::CommandLine::Option> opts{
    {'A', "accuracy", "RELERR", "", "", "",
     "Maximum relative error for least-squares fit (default: 1e-2)."},

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

    {'c', "check", "", "", "", "",
     "Generate comparative PACS color correction tables and exit."},

    {'d', "dust-model", "NAME", "", "", "",
     "Dust opacity model (default: OH5)."},

    {'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', "gas-to-dust", "RATIO", "", "", "",
     "Gas-to-dust ratio (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', "max-iterations", "NMAX", "", "", "",
     "Maximum least-squares iterations per pixel (default: 25)."},

    {'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', "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.)"},

    {'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('c')) {
    manticore::ccfCheck();
    exit(EXIT_SUCCESS);
  }

  manticore::process(cli);

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