#ifndef _MATERIAL_PARAMETERS_H_
#define _MATERIAL_PARAMETERS_H_

#include <cmath>
#include <utility>

/**
 * We indicate the test case we want to run
 */
enum TestCase{
  TEST_CASE_MOVE, ///< Drop movement
  TEST_CASE_SPLIT, ///< Drop splitting
  TEST_CASE_MERGE ///< Drop merging
};

/**
 * A material function that depends on the phase.
 *
 * In other words, if $m_1$ and $m_2$ are the values at the pure phases,
 * this is a function defined as:
 * \f[
 *   m( x ) = (m_1 - m_2 )\chi_{\{ \phi(x) >0 \}} + m_2
 * \f]
 */
struct material_function{
  public:
    /**
     * Explicit constructor.
     * @param _m1 :  value for $ x = 1 $
     * @param _m2 :  value for $x = -1$
     * @param _delta : interface thickness
     */
    material_function( const double _m1 = 1., const double _m2 = 1., const double _delta = 1. );
    /**
     * Constructor from pair.
     * @param params :  the pair of values for $x=1$ and $x=-1$ respectively
     * @param _delta : interface thickness
     */
    material_function( const std::pair<double, double> &params, const double _delta );
    /**
     * Derivative of the function or a smoothed version of it
     */
    double deriv( const double x ) const;
    /**
     * Maximal value of the derivative
     */
    double max_deriv() const;
    /**
     * The value
     */
    double operator()( const double x ) const;
    /**
     * <b>operator</b>()(<code>x</code>) + <b>operator</b>()(<code>y</code>)
     */
    double average( const double x, const double y ) const;
  private:
    const double avg, ///< The average between the two values
                 diff, ///< half the difference between the two values
                 delta; ///< The interface thickness
    /**
     * The function that actually computes the value
     */
    double value ( const double x ) const;
};

/**
 * This class implements the (modified) Ginzburg Landau potential, which is given by
 * \f[
 *    \mathcal{W}(x) = 
        \begin{cases}
          (1+x)^2, & x \le -1, \\
          \frac14 (1-x^2)^2, & |x| \leq 1, \\
          (1-x)^2, & x > 1.
        \end{cases}
 * \f]
 */
struct GLPotential{
  public:
    /**
     * Constructor. Just initalize the members
     */
    GLPotential();
    /**
     * Derivative: $\mathcal{W}'(x)$.
     */
    double operator()( const double x ) const;
  private:
    /**
     *  Maximum of the second derivative: $\sup_{x \in \mathbb{R}} |\mathcal{W}''(x)|$.
     */
    const double max_second_deriv;
  public:
    /**
     * The value of the maximum of the second derivative is known beforehands, and so
     * it is the value of the stabilization parameter. Thus, we make it a constant.
     */
    const double A;
};

/**
 * The Interface Energy density function, which is given by
 * \f[
     \gamma_{fs}(x) = \gamma \frac12 \cos \theta_s \sin(\frac{\pi x}2).
 * \f]
 * Notice that this function differs from the one used in the problem definition.
 * The relation is
 * \f[
 *    \gamma_{fs}(x) = \gamma \Theta_{fs}(x),
 * \f]
 * where $\gamma$ is the surface tension.
 */
struct InterfaceEnergy{
  public:
    /**
     * Initialize the member variables
     * @param gamma : The surface tension
     * @param theta_s : The static contact angle
     */
    InterfaceEnergy( const double gamma, const double theta_s );
    /**
     * Derivative: $\gamma_{fs}'(x)$.
     */
    double operator()( const double x ) const;
  private:
    const double coeff, ///< The coefficient that goes in front of the phase dependent function
                 max_second_deriv; ///<  $\sup_{x \in \mathbb{R}} |\gamma_{fs}''(x) |$.
  public:
    /**
     * The value of the maximum of the second derivative is known beforehand, and so
     * it is the value of the stabilization parameter. Thus, we make it a constant.
     */
    const double B;
};

/**
 * We bundle all the material parameters together. These are:
 * <ul>
 *  <li> $\varepsilon$, the permittivity of the fluids
 *  <li> $K$ the conductivity of the fluids
 *  <li> $M$, the permittivity of the fluids
 *  <li> $\rho$, the density of the fluids
 *  <li> $\beta$, the permittivity of the fluids
 *  <li> $\eta$, the permittivity of the fluids
 *  <li> $\lambda$, regularization parameter
 *  <li> $\gamma$, the surface tension
 *  <li> $\delta$, the interface thickness
 *  <li> $\alpha$, the surface mobility parameter
 *  <li> $\varepsilon^\star$, the conductivity of the plates
 *  <li> $\mathcal{W}$, the Ginzburg Landau potential
 *  <li> $\gamma_{fs}$, the interface energy density
 * </ul>
 * The two constructors are meant to get the information from either a bunch of numbers
 * or from pairs.
 */
struct Material_Parameters{
  public:
    Material_Parameters( const double epsilon1, const double epsilon2, const double epsilon3,
                         const double K1, const double K2, const double M1, const double M2,
                         const double rho1, const double rho2, const double beta1, const double beta2,
                         const double eta1, const double eta2, const double _lambda, const double _gamma,
                         const double theta_s, const double _delta, const double _alpha );
    Material_Parameters( std::pair<double, double> _epsilon, const double epsilon3,
                         std::pair<double, double> _K, 
                         std::pair<double, double> _M, std::pair<double, double> _rho,
                         std::pair<double, double> _beta, 
                         std::pair<double, double> _eta, const double _lambda, const double _gamma,
                         const double theta_s, const double _delta, const double _alpha );
    material_function epsilon, K, M, rho, beta, eta;
    const double lambda, gamma, delta, alpha, epsilon_plate;
    GLPotential W;
    InterfaceEnergy gamma_fs;
};

#endif //_MATERIAL_PARAMETERS_H_
