#ifndef TANH4_EX1_H //Usual macro guard to prevent multiple inclusion
#define TANH4_EX1_H

#include "CosmoInterface/cosmointerface.h"

//Include cosmointerface to have access to all of the library.

namespace TempLat
{
    /////////
    // Model name and number of fields
    /////////

    // In the following class, we define the defining parameters of your model:
    //number of fields of each species and the type of interactions.

    struct ModelPars : public TempLat::DefaultModelPars {
        static constexpr size_t NScalars = 1;
        static constexpr size_t NPotTerms = 1;
        static constexpr size_t NDim = 2;   
    };

  #define MODELNAME tanh4_Ex1
  // Here we define the name of the model. This should match the name of your file.

  template<class R>
  using Model = MakeModel(R, ModelPars);
  // In this line, we define an appropriate generic model, with the correct
  // number of fields, ready to be customized.
  // If you are curious about what this is doing, the macro is defined in
  // the "cosmointerface/abstractmodel.h" file.

  class MODELNAME : public Model<MODELNAME>
  //Declaration of our model. It inherits from the generic model defined above.
  {
 //...
private:

	double M, Lambda4, phii, omega, q;


  public:
  
      static constexpr int NPot = 4;

    MODELNAME(ParameterParser& parser, RunParameters<double>& runPar, std::shared_ptr<MemoryToolBox> toolBox): //Constructor of our model.
    Model<MODELNAME>(parser,runPar.getLatParams(), toolBox, runPar.dt, STRINGIFY(MODELLABEL)) //MODELLABEL is defined in the cmake.
    {
      /////////
      // Independent parameters of the model and initial homogeneous components of the fields
      /////////
      
        M = parser.get<double>("M");
        Lambda4 = parser.get<double>("Lambda4");

        
        fldS0 = parser.get<double, Ns>("initial_amplitudes");
        piS0 = parser.get<double, Ns>("initial_momenta");
        
    	phii = fldS0[0];
    	omega = sqrt(Lambda4 * pow<NPot-2>(phii) / pow<NPot>(M));

        
        /////////
        // Rescaling for program variables
        /////////
        
   		 fStar = phii;
   		 omegaStar = omega;
   		 alpha = 3.0 * (NPot - 2.0) / (NPot + 2.0);
   		 
   		 // We now need to specify the rescaling from physical units to program units.
        // This consists of the  time rescaling exponent alpha, the field rescaling fStar
        // and the velocity rescaling omegaStar.
        // See the paper for more information on how to fix them.
        

   		 setInitialPotentialAndMassesFromPotential();
   		 // Here we call this function to compute the value of the potential on the homogeneous
        // initial condition  (useful to set the initial Hubble rate). We also compute
        // in this function the masses from the second derivative of the potential
        // evaluated on the homogeneous initial conditions. If you want to do something else,
        // uncomment the next section and do whatever suits your needs.

        /*
          masses2S = {..., ...};
          setInitialPotentialFromPotential();
         */ 
   		
   	}
   		 
   /////////
   // Program potential (add as many functions as terms are in the potential)
   /////////  
   
    auto potentialTerms(Tag<0>) // Inflaton potential energy
    {
      return /*...*/ ;
    }

    
   /////////
   // Derivatives of the program potential with respect fields (add one function for each field)
   ///////// 

    auto potDeriv(Tag<0> f0) // Derivative with respect inflaton
    {
        return  /*...*/ ;
    }
    
   /////////
   //  Second derivatives of the program potential with respect fields (add one function for each field)
   /////////

    auto potDeriv2(Tag<0> f0) // Second derivative with respect inflaton
    {
        return  /*...*/ ;
    }



    };
}

#endif //TANH4_EX1_H