#ifndef __JMATH__JMATHLIB3D__
#define __JMATH__JMATHLIB3D__

#include "JMath/JMathlib.hh"

/**
 * \file
 * Functional algebra in 3D.
 *
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  /**
   * Make 3D function of x from 1D function.
   */
  template<class JF1_t>
  struct JMake3X :
    public JMathlib< JMake3X<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JMake3X<JF1_t> >::operator();


    /**
     * Default constructor.
     */
    JMake3X()
    {}


    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JMake3X(const JF1_t& f1) :
      JF1_t(f1)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JMake3X(const Args& ...args) :
      JF1_t(args...)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y, const double z) const
    {
      return static_cast<const JF1_t&>(*this).getValue(x);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               gradient
     */
    const JMake3X& getGradient(const double x, const double y, const double z) const
    {
      static JMake3X gradient;

      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(x);

      return gradient;
    }
  };


  /**
   * Make 3D function of y from 1D function.
   */
  template<class JF1_t>
  struct JMake3Y :
    public JMathlib< JMake3Y<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JMake3Y<JF1_t> >::operator();


    /**
     * Default constructor.
     */
    JMake3Y()
    {}


    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JMake3Y(const JF1_t& f1) :
      JF1_t(f1)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JMake3Y(const Args& ...args) :
      JF1_t(args...)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y, const double z) const
    {
      return static_cast<const JF1_t&>(*this).getValue(y);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               gradient
     */
    const JMake3Y& getGradient(const double x, const double y, const double z) const
    {
      static JMake3Y gradient;

      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(y);

      return gradient;
    }
  };


  /**
   * Make 3D function of z from 1D function.
   */
  template<class JF1_t>
  struct JMake3Z :
    public JMathlib< JMake3Z<JF1_t> >,
    public JF1_t
  {
    using JMathlib< JMake3Z<JF1_t> >::operator();


    /**
     * Default constructor.
     */
    JMake3Z()
    {}


    /**
     * Constructor.
     *
     * \param  f1            function
     */
    JMake3Z(const JF1_t& f1) :
      JF1_t(f1)
    {}


    /**
     * Constructor.
     *
     * \param  args          list of values
     */
    template<class ...Args>
    JMake3Z(const Args& ...args) :
      JF1_t(args...)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y, const double z) const
    {
      return static_cast<const JF1_t&>(*this).getValue(z);
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               gradient
     */
    const JMake3Z& getGradient(const double x, const double y, const double z) const
    {
      static JMake3Z gradient;

      static_cast<JF1_t&>(gradient) = static_cast<const JF1_t&>(*this).getGradient(z);

      return gradient;
    }
  };


  template<int ID_t, size_t N> using JPolynome3X = JMake3X< JPolynome<ID_t, N> >;                    //!< 3D polynomial function of x.
  template<int ID_t, size_t N> using JPolynome3Y = JMake3Y< JPolynome<ID_t, N> >;                    //!< 3D polynomial function of y.
  template<int ID_t, size_t N> using JPolynome3Z = JMake3Z< JPolynome<ID_t, N> >;                    //!< 3D polynomial function of z.
  template<int ID_t, bool normalise = false> using JGauss3X = JMake3X< JGauss<ID_t, normalise> >;    //!< 3D Gauss function of x.
  template<int ID_t, bool normalise = false> using JGauss3Y = JMake3Y< JGauss<ID_t, normalise> >;    //!< 3D Gauss function of y.
  template<int ID_t, bool normalise = false> using JGauss3Z = JMake3Z< JGauss<ID_t, normalise> >;    //!< 3D Gauss function of z.
  template<int ID_t, class JF1_t = _vF> using JPow3X = JMake3X< JPow<ID_t, JF1_t> >;                 //!< 3D power of function of x.
  template<int ID_t, class JF1_t = _vF> using JPow3Y = JMake3Y< JPow<ID_t, JF1_t> >;                 //!< 3D power of function of y.
  template<int ID_t, class JF1_t = _vF> using JPow3Z = JMake3Z< JPow<ID_t, JF1_t> >;                 //!< 3D power of function of z.
  template<int N>  using JXn3X = JMake3X< JXn<N> >;                                                  //!< 3D fixed power of x.
  template<int N>  using JXn3Y = JMake3Y< JXn<N> >;                                                  //!< 3D fixed power of y.
  template<int N>  using JXn3Z = JMake3Z< JXn<N> >;                                                  //!< 3D fixed power of z.
  template<class JF1_t = _vF> using JSqrt3X = JMake3X< JSqrt<JF1_t> >;                               //!< 3D square root of function of x.
  template<class JF1_t = _vF> using JSqrt3Y = JMake3Y< JSqrt<JF1_t> >;                               //!< 3D square root of function of y.
  template<class JF1_t = _vF> using JSqrt3Z = JMake3Z< JSqrt<JF1_t> >;                               //!< 3D square root of function of z.
  template<class JF1_t> using JSin3X = JMake3X< JSin<JF1_t> >;                                       //!< 3D sine of function of x.
  template<class JF1_t> using JSin3Y = JMake3Y< JSin<JF1_t> >;                                       //!< 3D sine of function of y.
  template<class JF1_t> using JSin3Z = JMake3Z< JSin<JF1_t> >;                                       //!< 3D sine of function of z.
  template<class JF1_t> using JCos3X = JMake3X< JCos<JF1_t> >;                                       //!< 3D cosine of function of x.
  template<class JF1_t> using JCos3Y = JMake3Y< JCos<JF1_t> >;                                       //!< 3D cosine of function of y.
  template<class JF1_t> using JCos3Z = JMake3Z< JCos<JF1_t> >;                                       //!< 3D cosine of function of z.
  template<class JF1_t> using JExp3X = JMake3X< JExp<JF1_t> >;                                       //!< 3D exponent of function of x.
  template<class JF1_t> using JExp3Y = JMake3Y< JExp<JF1_t> >;                                       //!< 3D exponent of function of y.
  template<class JF1_t> using JExp3Z = JMake3Z< JExp<JF1_t> >;                                       //!< 3D exponent of function of z.
  template<class JF1_t> using JLog3X = JMake3X< JLog<JF1_t> >;                                       //!< 3D logarithm of function of x.
  template<class JF1_t> using JLog3Y = JMake3Y< JLog<JF1_t> >;                                       //!< 3D logarithm of function of y.
  template<class JF1_t> using JLog3Z = JMake3Z< JLog<JF1_t> >;                                       //!< 3D logarithm of function of z.

  
  /**
   * 3D correlated Gauss function.
   */
  template<int ID_t, bool normalise = false>
  struct JGauss3D :
    public JMathlib< JGauss3D<ID_t, normalise> >,
    public JGauss<ID_t, normalise>
  {
    using JMathlib< JGauss3D<ID_t, normalise> >::operator();
    

    /**
     * Default constructor.
     */
    JGauss3D()
    {}


    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     */
    JGauss3D(const double mean,
             const double sigma) :
      JGauss<ID_t, normalise>(mean, sigma)
    {}


    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               function value
     */
    double getValue(const double x, const double y, const double z) const
    {
      typedef JGauss<ID_t, normalise>  JGauss_t;

      return (static_cast<const JGauss_t&>(*this).getValue(x) *
              static_cast<const JGauss_t&>(*this).getValue(y) *
              static_cast<const JGauss_t&>(*this).getValue(z));
    }


    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \param  y             abscissa value
     * \param  z             abscissa value
     * \return               gradient
     */
    const JGauss3D& getGradient(const double x, const double y, const double z) const
    {
      static JGauss3D gradient;

      typedef JGauss<ID_t, normalise>  JGauss_t;

      static_cast<JGauss_t&>(gradient)  = JGauss_t(static_cast<const JGauss_t&>(*this).getGradient(x)).mul(static_cast<const JGauss_t&>(*this).getValue(y) * static_cast<const JGauss_t&>(*this).getValue(z));
      static_cast<JGauss_t&>(gradient) += JGauss_t(static_cast<const JGauss_t&>(*this).getGradient(y)).mul(static_cast<const JGauss_t&>(*this).getValue(x) * static_cast<const JGauss_t&>(*this).getValue(z));
      static_cast<JGauss_t&>(gradient) += JGauss_t(static_cast<const JGauss_t&>(*this).getGradient(z)).mul(static_cast<const JGauss_t&>(*this).getValue(x) * static_cast<const JGauss_t&>(*this).getValue(y));

      return gradient;
    }
  };
}

#endif