#ifndef DUST_H
#define DUST_H

#include <gsl/gsl_spline.h>

#include "mutils/CommandLine.h"
#include "mutils/util.h"

namespace manticore {

namespace dust {
  extern void printModels(std::vector<mutils::CommandLine::Option> &opts,
                          char shortName);
}


/**
  \brief Dust model.

  This class supports both arbitrary, table-based opacity models as well as
  fixed power laws. The underlying model is set at construction but can be
  changed during use; cf. setModel(). Table-based models are evaluated via
  cubic spline interpolation of the logarithm of the opacity.

  \note
  Once initialized, a model may be safely shared among multiple threads.
*/
class Dust {
public:
  /// Default constructor.
  /// \param[in] model Dust model name.
  /// \param[in] rho   Gas-to-dust ratio.
  Dust(const std::string &model = "OH5", double rho = 100.0)
  { setModel(model, rho); }

  /// Destructor.
  ~Dust() {
    if (spline_) { gsl_spline_free(spline_); }
  }

  /// Fixed power-law constructor.
  /// \param[in] nu0     Fiducial frequency (Hz).
  /// \param[in] kappa0  Fiducial opacity (cm^2/g @ g/d = 100).
  /// \param[in] beta   %Opacity scaling.
  Dust(double nu0, double kappa0, double beta)
  { setModel(nu0, kappa0, beta); }

  /// Set dust properties model (table-based).
  /// \param[in] model Dust model name.
  /// \param[in] rho   Gas-to-dust ratio.
  void setModel(const std::string &model, double rho = 100.0);

  /// Set dust properties model (fixed power law).
  /// \param[in] nu0     Fiducial frequency (Hz).
  /// \param[in] kappa0  Fiducial opacity (cm^2/g @ g/d = 100).
  /// \param[in] beta   %Opacity scaling.
  void setModel(double nu0, double kappa0, double beta);

  /// Gas/dust extinction opacity (cm^2/g).
  /// \param[in] nu  Observation frequency (Hz).
  /// \param[in] acc Interpolation lookup accelerator.
  double kappa(double nu, gsl_interp_accel *acc = nullptr) const;

  /// Current model name.
  const std::string &name() const noexcept { return name_; }

  /// Current model gas-to-dust ratio.
  double rho() const noexcept { return rho_; }

protected:
  /// Fixed power-law opacity parameters.
  struct Opacity { double inu0, kappa0, beta; };

  /// Current model name.
  std::string name_;

  /// Current model gas-to-dust ratio.
  double rho_;

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

  /// Table interpolator.
  gsl_spline *spline_ = nullptr;

  /// Whether to interpolate in log space.
  bool ln_ = false;
};

}  // namespace manticore

#endif  // DUST_H