/**
  \file

  \brief Dust model specifications.

  \who \kpr
*/
#include <cmath>

#include "Dust.h"
#include "manticore.h"

namespace mu = mutils;


namespace manticore {

/// Dust models.
namespace dust {

/**
  \brief OH5 model.

  Ossenkopf & Henning (1994) thinly ice-mantled coagulated dust, augmented by
  Pollack et al. (1994) opacities at < 1.25 microns (as described in Dunham et
  al. 2010).
*/
const manticore::tableType modelOH5 =
{ // {nu0, kappa0}
  {8.321743e+09, 2.286790e-05},
  {1.248261e+10, 4.744500e-05},
  {1.872392e+10, 9.843610e-05},
  {2.304483e+10, 1.430450e-04},
  {2.723480e+10, 1.932250e-04},
  {4.048416e+10, 3.944150e-04},
  {5.991655e+10, 7.987750e-04},
  {8.559507e+10, 1.517920e-03},
  {8.942755e+10, 1.642440e-03},
  {1.097373e+11, 2.373980e-03},
  {1.152241e+11, 2.591900e-03},
  {1.331476e+11, 3.362340e-03},
  {1.497914e+11, 4.156390e-03},
  {1.997218e+11, 6.975990e-03},
  {2.304483e+11, 9.025520e-03},
  {2.995828e+11, 1.369980e-02},
  {4.279754e+11, 2.569960e-02},
  {5.991655e+11, 5.039920e-02},
  {8.559507e+11, 1.009980e-01},
  {1.325589e+12, 2.169970e-01},
  {1.621122e+12, 3.076000e-01},
  {1.890111e+12, 4.066510e-01},
  {2.204437e+12, 5.256380e-01},
  {2.569312e+12, 6.786530e-01},
  {2.995828e+12, 8.649890e-01},
  {3.492859e+12, 1.180770e+00},
  {4.072078e+12, 1.581190e+00},
  {4.747743e+12, 2.109970e+00},
  {5.535526e+12, 2.798060e+00},
  {6.453742e+12, 3.866420e+00},
  {6.968660e+12, 4.210170e+00},
  {7.525308e+12, 4.269800e+00},
/* lambda < 40 micron:
  {8.125382e+12, 4.341700e+00},
  {8.772556e+12, 4.506770e+00},
  {1.022816e+13, 5.454030e+00},
  {1.192606e+13, 7.666310e+00},
  {1.390174e+13, 1.094070e+01},
  {1.501665e+13, 1.347080e+01},
  {1.621122e+13, 1.579020e+01},
  {1.750920e+13, 1.510520e+01},
  {1.890111e+13, 1.336780e+01},
  {2.040751e+13, 1.354000e+01},
  {2.204437e+13, 1.603100e+01},
  {2.379532e+13, 1.965150e+01},
  {2.569312e+13, 2.513110e+01},
  {2.670067e+13, 2.548350e+01},
  {2.773904e+13, 2.480020e+01},
  {2.883368e+13, 2.521440e+01},
  {2.995828e+13, 2.579990e+01},
  {3.053850e+13, 2.679960e+01},
  {3.113194e+13, 2.207110e+01},
  {3.234883e+13, 1.721290e+01},
  {3.361565e+13, 1.431590e+01},
  {3.492859e+13, 1.396440e+01},
  {3.771657e+13, 1.331150e+01},
  {4.072631e+13, 1.268920e+01},
  {4.397222e+13, 1.209670e+01},
  {4.747743e+13, 1.153170e+01},
  {5.127208e+13, 1.369520e+01},
  {5.535526e+13, 1.290380e+01},
  {5.977308e+13, 1.419280e+01},
  {6.453742e+13, 1.628710e+01},
  {6.968660e+13, 1.760470e+01},
  {7.525308e+13, 1.929400e+01},
  {7.819966e+13, 2.029280e+01},
  {8.125382e+13, 2.133790e+01},
  {8.443693e+13, 2.322900e+01},
  {8.775124e+13, 2.533260e+01},
  {9.116944e+13, 2.974520e+01},
  {9.474451e+13, 4.313220e+01},
  {9.844981e+13, 6.182070e+01},
  {1.022816e+14, 6.525970e+01},
  {1.063105e+14, 3.606000e+01},
  {1.104656e+14, 3.394000e+01},
  {1.192606e+14, 3.894270e+01},
  {1.390819e+14, 5.063830e+01},
  {1.621122e+14, 6.449570e+01},
  {1.890111e+14, 7.974990e+01},
  {2.204437e+14, 9.807170e+01},
  {2.569312e+14, 1.143450e+02},
  {2.995828e+14, 1.309990e+02},
  {3.153503e+14, 1.370850e+02},
  {3.328697e+14, 1.438060e+02},
  {3.524503e+14, 1.512710e+02},
  {3.744784e+14, 1.596120e+02},
  {3.994437e+14, 1.689970e+02},
  {4.279754e+14, 1.796410e+02},
  {4.608965e+14, 1.918230e+02},
  {4.993046e+14, 2.059090e+02},
  {5.446959e+14, 2.223980e+02},
  {5.991655e+14, 2.419780e+02},
  {6.657394e+14, 2.656360e+02},
  {7.489569e+14, 2.948310e+02},
  {8.559507e+14, 3.318290e+02},
  {9.986092e+14, 3.803510e+02},
  {1.198331e+15, 4.469770e+02},
  {1.497914e+15, 5.446050e+02},
  {1.997218e+15, 7.025750e+02},
  {2.995828e+15, 1.005980e+03},
  {5.991655e+15, 1.858230e+03},
*/
};


// Register new models here as needed.
std::map<std::string, manticore::tableType>
modelMap =
{
  {"OH5", dust::modelOH5},
};


/**
  \brief Add available dust model names to options summary.
  \param[in] opts Command line options.
  \param[in] name Short name of dust models option.
*/
void printModels(std::vector<mu::CommandLine::Option> &opts, char name)
{
  // Find dust model option.
  for (auto &entry: opts) {
    if (entry.shortName == name) {
      entry.helpText += "\n(Available models:";
      for (auto &model: modelMap) { entry.helpText += " " + model.first; }
      entry.helpText += ")";
    }
  }
}

}  // namespace dust



void Dust::setModel(const std::string &model, double rho)
{
  const char *const fn = "Dust::setModel";
  auto mp = dust::modelMap.find(model);
  if (mp == dust::modelMap.end()) {
    mp = dust::modelMap.begin();
    MU_ERROR(fn, "Unknown dust model '%s', using '%s'", model, mp->first);
  }

  name_  = mp->first;
  rho_   = rho;

  auto table = mp->second;
  if (rho != 100.0) {
    double scale = 100.0/rho;
    for (auto &entry: table) { entry.second *= scale; }
  }

  // Release old model.
  if (spline_) {
    gsl_spline_free(spline_);
    spline_ = nullptr;
  }

  // Interpolate logarithm of ordinates.
  spline_ = initSpline(table, true);

  // Invalidate fixed power law.
  opacity_ = Opacity{0.0, 0.0, 0.0};
}


void Dust::setModel(double lamb0, double kappa0, double beta)
{
  char name[32];
  snprintf(name, sizeof(name), "beta=%.6g", beta);
  name_ = name;

  opacity_ = Opacity{lamb0/(1e4*c_light), kappa0, beta};

  // Ratio not meaningful here.
  rho_ = 0.0;

  // Free old model.
  if (spline_) {
    gsl_spline_free(spline_);
    spline_ = nullptr;
  }
}


/**
  The opacity value is interpolated from tabulated values assuming a fixed
  power law dependence between adjacent points. If \a nu lies woutside the
  range of the table, the corresponding end value is extrapolated.

  \return Dust+gas extinction opacity at the given wavelength.
*/
double Dust::kappa(double nu) const
{
  if (spline_) {
    return exp(gsl_spline_eval(spline_, nu, nullptr));
  } else {
    // Fixed power law.
    return opacity_.kappa0 * pow(nu*opacity_.inu0, opacity_.beta);
  }
}

}  // namespace manticore