#include "../include/utilities.h"
#include "../include/ParametersManager.h"

extern "C" {

//------------------------------------------------------------------------------------------------------ 
//  Model description header  for the problem Lac
//  Edit this file to define your model parameters 
//------------------------------------------------------------------------------------------------------ 
std::string paramsFile = "Lake_params.json";
//------------------------------------------------------------------------------------------------------ 
// Dynamics definition  
// Please fill the code of the below function without changing it's signature  
//------------------------------------------------------------------------------------------------------ 

double b = 0.8;
double q = 8;
double mm = 1;
double r = 1;
double mq = pow(mm,q);
double Pmax = 1.2;
double Lmin = 0.1;
double Lmax = 1;
double umax = 0.1;

void dynamics(const double *x, const double *u, double *image)
{ 
// programm here your dynamics function: 
// x :     the input  state vector
// u :     the input control vector 
// image : the resulting state vector 
double pq = x[1] > 0 ? pow(x[1], q) : 0.0;
image[0]= u[0];
image[1]= - b * x[1] + x[0] + r * pq / (mq + pq);
} 
//------------------------------------------------------------------------------------------------------ 
// Jacobian matrix of the dynamics. It will be used to estimate locally the Lipschitz constant  
// Please fill the code of the below function without changing it's signature  
//------------------------------------------------------------------------------------------------------ 
void jacobian(const double *x, const double *u , double ** jacob){
// programm here the jacobian matrix of your dynamics function: 
// x :     the input  state vector
// u :     the input control vector 
// jacob : the resulting matrix 
double pq = x[1] > 0 ? pow(x[1], q) : 0.0;
double pq1 = x[1] > 0 ? pow(x[1], q - 1.0) : 0.0;
jacob[0][0]=0.0;
jacob[0][1]=0.0;
jacob[1][0]=1.0;
jacob[1][1]=-b + mq * pq1 / ((mq + pq) * (mq + pq));
} 
//------------------------------------------------------------------------------------------------------ 
// Local bound of the dynamics. Each element of the result vector is defined as max(over u) F_i(x,u)  
// Please fill the code of the below function without changing it's signature  
//------------------------------------------------------------------------------------------------------ 
void localDynBounds(const double *x, double * bound){
// programm here the bound of your dynamics function, maximized over the control: 
// x :     the input  state vector
// bound : the resulting vector 
bound[0]=umax;
bound[1]=b * Pmax + Lmax + 1.0;
} 


//------------------------------------------------------------------------------------------------------ 
// Function defining the constraints set K. It should return a finite double value for any point  
// inside K, and the infinity for any point outside K  
// IMPORTANT : VIABLAB defines infinity as constant PLUS_INF. Please use it in this function  
// IMPORTANT : for viability problems with epigraphic property, this function is used as 
// initial state of the value function 
// Please fill the code of teh below function without changing it's signature  
//------------------------------------------------------------------------------------------------------ 


 double constraintsX( const double *x )
{
// programm here your constraints function: 
// x : the input  state vector
// The function returns PLUS_INF for points outside the contrants set
// and a double value for points inside the constraints set


bool outOfK = false;
outOfK =  (x[0]<Lmin) ||  (x[0]>Lmax) || (x[1]>Pmax) || (x[1]<0);
 
 return outOfK ? PLUS_INF : 1.0;
}


}