CDPP / AMDA_Kernel

Blame view

src/ExternLib/Deriv/Deriv.hh 3.98 KB
Edit Raw Normal View History
fbe3c2bb   Benjamin Renard   First commit 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
 
/*
 * 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_ */