/*
 * Deriv.hh
 *
 *  Created on: Dec 3, 2012
 *      Author: f.casimir
 */

#ifndef Deriv_HH_
#define Deriv_HH_

#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"

namespace AMDA {
namespace Parameters {

/**
 * @class Deriv
 * @brief It is responsible to compute the deriv data.
 * @details This class implement the interface Operation.
 * @code
 * 		3 point Lagrangian interpolation (points {x0,y0; x1,y1; x2,y2}):
 * 		df/dx = y0*(2x-x1-x2)/(x01*x02) + y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21),
 * 		where
 * 		x01= x0-x1; x02=x0-x2, x12=x1-x2...
 * @endcode
 */
template<class TParamData>
class Deriv : public Operation {
public:
	/**
	 * @brief Constructor.
	 * @details Create the ParamData type of the input ParamData.
	 */
  Deriv(Process& pProcess, TParamData& paramInput)
	: Operation(pProcess),
	  _paramInput(paramInput),
	  _paramOutput(new TParamData()),
	  _index(0),
	  _lastDataWrite(false) {
	  _paramDataOutput=_paramOutput;
  }

	/**
	 * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
	 */
  void  write(ParamDataIndexInfo &pParamDataIndexInfo) {

    typename TParamData::ElementType y0;
    typename TParamData::ElementType y1;
    typename TParamData::ElementType y2;
    typename TParamData::ElementType res;
    double x0 = 0;
    double x1 = 0;
    double x2 = 0;
    double x02 = 0;
    double x01 = 0;
    double x10 = 0;
    double x12 = 0;
    double x20 = 0;
    double x21 = 0;
    double _2x = 0;
    double x = 0;
    //if first compute
    if(_index == 0 ) {
    	if ((pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess) > 2) {
			y0 = _paramInput.getDataList()[_index];
			y1 = _paramInput.getDataList()[_index + 1];
			y2 = _paramInput.getDataList()[_index + 2];
			x = x0 =_paramInput.getTime(_index) ;
			x1 = _paramInput.getTime(_index + 1) ;
			x2 = _paramInput.getTime(_index + 2) ;
			x02 = x0 -x2;
			x01 = x0 -x1;
			x10 = x1 -x0;
			x12 = x1 -x2;
			x20 =x2 - x0 ;
			x21 = x2 -x1;
			_2x = 2*x;
			res = y0*(_2x-x1-x2)/(x01*x02) + y1*(_2x-x0-x2)/(x10*x12)+y2*(_2x-x0-x1)/(x20*x21);
		    _paramOutput->getDataList().push_back(res);
			_paramOutput->pushTime(x);
			++_index;
		}
    }
    //if last compute
     if( pParamDataIndexInfo._nbDataToProcess == 0 && !_lastDataWrite) {

		y0 = _paramInput.getDataList()[_index - 2];
		y1 = _paramInput.getDataList()[_index - 1];
		y2 = _paramInput.getDataList()[_index ];
		x0 =_paramInput.getTime(_index - 2) ;
		x1 = _paramInput.getTime(_index - 1) ;
		x = x2 = _paramInput.getTime(_index) ;
		x02 = x0 -x2;
		x01 = x0 -x1;
		x10 = x1 -x0;
		x12 = x1 -x2;
		x20 =x2 - x0 ;
		x21 = x2 -x1;
		_2x = 2*x;
		res = y0*(_2x-x1-x2)/(x01*x02) + y1*(_2x-x0-x2)/(x10*x12)+y2*(_2x-x0-x1)/(x20*x21);
	    _paramOutput->getDataList().push_back(res);
		_paramOutput->pushTime(x);
		_lastDataWrite = true;

    }
    //classic compute
	while (_index< pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess - 1) {
			//df/dt = y0*(2x-x1-x2)/(x01*x02) + y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21),

				y0 = _paramInput.getDataList()[_index - 1];
				y1 = _paramInput.getDataList()[_index ];
				y2 = _paramInput.getDataList()[_index + 1];
				x0 =_paramInput.getTime(_index - 1) ;
				x = x1 = _paramInput.getTime(_index) ;
				x2 = _paramInput.getTime(_index + 1) ;
				x02 = x0 -x2;
				x01 = x0 -x1;
				x10 = x1 -x0;
				x12 = x1 -x2;
				x20 =x2 - x0 ;
				x21 = x2 -x1;
				_2x = 2*x;
				res = y0*(_2x-x1-x2)/(x01*x02) + y1*(_2x-x0-x2)/(x10*x12)+y2*(_2x-x0-x1)/(x20*x21);
			    _paramOutput->getDataList().push_back(res);
				_paramOutput->pushTime(x);
				++_index;
		}

  }

private:
	/**<
	 * @brief It is the channel of data derived
	 */
  TParamData &_paramInput;
	/**<
	 * @brief It is the channel of the data derived
	 */
  TParamData *_paramOutput;

  /**<
   * index current of compute
   */
  unsigned int _index ;

  /**<
   * true when last data is write.
   */
  bool _lastDataWrite;
};

} /* namespace Parameters */
} /* namespace AMDA */
#endif /* Deriv_HH_ */