#include "AnT.hpp"

#define a       parameters[0]
#define gamma   parameters[1]
#define XL      currentState[cellIndex-1][0]
#define X       currentState[cellIndex][0]
#define XR      currentState[cellIndex+1][0]
#define X_0     currentState[0][0]
#define X_N     currentState[numberOfCells-1][0]
#define XrL     referenceState[cellIndex-1][0]
#define Xr      referenceState[cellIndex][0]
#define XrR     referenceState[cellIndex+1][0]
#define Xr_0    referenceState[0][0]
#define Xr_N    referenceState[numberOfCells-1][0]

bool cll_linearized (const CellularState& currentState,
		     const CellularState& referenceState,
		     const Array<real_t>& parameters,
		     int cellIndex,
		     StateCell& rhs)
{
  int numberOfCells = currentState.numberOfCells;

  if (cellIndex == 0) {
      rhs[0] = (1-gamma)*(-2 * a * Xr * X)
	  +(gamma/2)*((-2 * a * Xr_N * X_N)+(-2 * a * XrR * XR));
  } else if (cellIndex == (numberOfCells-1)) {
      rhs[0] = (1-gamma)*(-2 * a * Xr * X)
	  +(gamma/2)*((-2 * a * XrL * XL)+(-2 * a * Xr_0 * X_0));
  }
  else {
      rhs[0] = (1-gamma)*(-2 * a * Xr * X)
	  +(gamma/2)*((-2 * a * XrL * XL)+(-2 * a * XrR * XR));
  }

  return true;
}



bool cll (const CellularState& currentState,
	  const Array<real_t>& parameters,
	  int cellIndex,
	  StateCell& rhs)
{
    int numberOfCells = currentState.numberOfCells;

    if (cellIndex == 0) {
	rhs[0] = (1-gamma)*(1-a*X*X)+(gamma/2)*((1-a*X_N*X_N)+(1-a*XR*XR));
    } else if (cellIndex == (numberOfCells-1)) {
	rhs[0] = (1-gamma)*(1-a*X*X)+(gamma/2)*((1-a*XL*XL)+(1-a*X_0*X_0));
    }
    else {
	rhs[0] = (1-gamma)*(1-a*X*X)+(gamma/2)*((1-a*XL*XL)+(1-a*XR*XR));
    }
    return true;
}

#undef alpha
#undef gamma
#undef XL    
#undef X     
#undef XR    
#undef X_0   
#undef X_N   
#undef XrL 
#undef Xr  
#undef XrR 
#undef Xr_0
#undef Xr_N

extern "C" 
{  
  void connectSystem ()
  {
    CML_Proxy::systemFunction = cll;
    CML_LinearizedProxy::systemFunction = cll_linearized;
  } 
}