//
// fermi.hpp
//
// LGPL Version 2.1 HEADER START
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301  USA
//
// LGPL Version 2.1 HEADER END
//

//
// Copyright (c) 2021, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Yaser Afshar
//

#ifndef FERMI_HPP
#define FERMI_HPP

#include <cmath>

#include "special.hpp"

/*! Fermi-type switching function */
class FERMI {
 public:
  /*!
   * \brief Construct a new FERMI object
   *
   * \param cutoff Cutoff distance
   * \param skin Width of the switching region.
   * \param shift Shift \f$ \delta \f$
   *
   */
  FERMI(double const cutoff, double const skin, double const shift);

  /*!
   * \brief Convert units of the parameters
   *
   * \param convert_length_factor length unit conversion factor
   * \param convert_energy_factor energy unit conversion factor
   */
  inline void ConvertUnit(double const convert_length_factor,
                          double const convert_energy_factor);

  /*!
   * \brief The Fermi switching function
   *
   * \f$f(r)=\left[1+\exp{\frac{r-r_c+\delta}{s}}\right]^{-1},\f$
   * where \f$r_c\f$ is the cutoff distance, \f$\delta\f$ is the shift and
   * \f$s\f$ is the skin (controls the sharpness of the transition).
   *
   * \param r distance
   * \param dfdr derivative of the Fermi switching function
   * \return double Fermi switching function value at r
   */
  inline double FermiSwitchingFunction(double const r, double &dfdr) const;

  /*!
   * \brief The CFermi switching function
   *
   * \f$c(r)=1-\left[1+\exp{\frac{r-r_c-\delta}{s}}\right]^{-1},\f$
   * where \f$r_c\f$ is the cutoff distance, \f$\delta\f$ is the shift and
   * \f$s\f$ is the skin (controls the sharpness of the transition).
   *
   * \param r distance
   * \param dfdr derivative of the CFermi switching function
   * \return double CFermi switching function value at r
   */
  inline double CFermiSwitchingFunction(double const r, double &dfdr) const;

 public:
  /*! Cutoff distance */
  double cutoff_;
  /*! Width of the switching region */
  double skin_;
  /*! Shift \f$ \delta \f$ */
  double shift_;
  /*! Skin inverse controls the sharpness of the transition. */
  double skin_inverse_;
};

inline FERMI::FERMI(double const cutoff, double const skin, double const shift)
    : cutoff_(cutoff), skin_(skin), shift_(shift) {
  skin_inverse_ = 1.0 / skin_;
}

inline void FERMI::ConvertUnit(double const convert_length_factor,
                               double const /* convert_energy_factor */) {
  if (special::IsNotOne(convert_length_factor)) {
    cutoff_ *= convert_length_factor;
    skin_ *= convert_length_factor;
    shift_ *= convert_length_factor;
    skin_inverse_ = 1.0 / skin_;
  }
}

inline double FERMI::FermiSwitchingFunction(double const r,
                                            double &dfdr) const {
  double const rcutoff = r - cutoff_ + shift_;
  double const rcutskin = rcutoff * skin_inverse_;
  double const fexp = std::exp(rcutskin);
  double const fexp1 = fexp + 1.0;
  double const denom = skin_ * fexp1 * fexp1;
  dfdr = -fexp / denom;
  double const f = 1.0 / fexp1;
  return f;
}

inline double FERMI::CFermiSwitchingFunction(double const r,
                                             double &dfdr) const {
  double const rcutoff = r - cutoff_ - shift_;
  double const rcutskin = rcutoff * skin_inverse_;
  double const fexp = std::exp(rcutskin);
  double const fexp1 = fexp + 1.0;
  double const denom = skin_ * fexp1 * fexp1;
  dfdr = fexp / denom;
  double const f = fexp / fexp1;
  return f;
}

#endif  // FERMI