CDPP / AMDA_Kernel

Blame view

src/ExternLib/Tsyganenko96/Tsyganenko96.hh 7.42 KB
Edit Raw Normal View History
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
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
 
/*
 * Tsyganenko96.hh
 *
 *  Created on: Nov 06, 2017
 *      Author: benjamin
 */

#ifndef TSYGANENKO96_HH_
#define TSYGANENKO96_HH_

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

#include <iterator>

namespace AMDA {
namespace Parameters {
namespace Tsyganenko96 {

/**
 * @class Tsyganenko96
 * @brief 
 * @details This class implement the interface Operation.
 */

#define AVERAGE_TIME 1200 // (seconds)
#define MAX_GAP_SIZE 3600 // (seconds)
3ab32f47   Benjamin Renard   Imf default value... 
31
32
33
 
#define DEFAULT_IMF_GSM_X 0.
#define DEFAULT_IMF_GSM_Y 2.
#define DEFAULT_IMF_GSM_Z -3.
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
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
 
#define DEFAULT_DST -10.
#define DEFAULT_PSW 3.

class Tsyganenko96Base : public Operation {
public:
	Tsyganenko96Base(Process& pProcess, ParamData& paramImfInput, ParamData& paramPswInput, ParamData& paramDstInput) 
	: Operation(pProcess),
	_paramImfInput(dynamic_cast<ParamDataSpec<std::vector<float> >&>(paramImfInput)),
	_paramPswInput(dynamic_cast<ParamDataSpec<float>&>(paramPswInput)),
	_paramDstInput(dynamic_cast<ParamDataSpec<float>&>(paramDstInput)) {
	}

	virtual ~Tsyganenko96Base() {
	}

	void pushImfData(ParamDataIndexInfo &pParamDataIndexInfo) {
		for (unsigned int _index = pParamDataIndexInfo._startIndex ;
			_index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
			++_index)
		{
			double time = _paramImfInput.getTime(_index);
			std::vector<float> inputElt = _paramImfInput.get(_index);

			_b_x_gse.push_back(std::pair<double,float>(time,inputElt[0]));
			_b_y_gse.push_back(std::pair<double,float>(time,inputElt[1]));
			_b_z_gse.push_back(std::pair<double,float>(time,inputElt[2]));
		}
	}

	void pushPswData(ParamDataIndexInfo &pParamDataIndexInfo) {
		for (unsigned int _index = pParamDataIndexInfo._startIndex ;
			_index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
			++_index)
		{
			double time = _paramPswInput.getTime(_index);
			float inputElt = _paramPswInput.get(_index);

			_psw.push_back(std::pair<double,float>(time,inputElt));
		}
	}

	void pushDstData(ParamDataIndexInfo &pParamDataIndexInfo) {
		for (unsigned int _index = pParamDataIndexInfo._startIndex ;
			_index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
			++_index)
		{
			double time = _paramDstInput.getTime(_index);
			float inputElt = _paramDstInput.get(_index);

			_dst.push_back(std::pair<double,float>(time,inputElt));
		}
	}
3ab32f47   Benjamin Renard   Imf default value... 
87
 
	float getValue(std::vector<std::pair<double,float> >& input, double time) {
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
 
		double min_t = time - AVERAGE_TIME/2.;
		double max_t = time + AVERAGE_TIME/2.;
		std::vector<std::pair<double,float> > values_for_mean;
		std::pair<double,float> prev_value(NAN,NAN);
		std::pair<double,float> next_value(NAN,NAN);
		for (std::vector<std::pair<double,float> >::iterator it = input.begin(); it != input.end(); ++it) {
			if (isNAN(it->second))
				continue;
			else if (it->first > max_t) {
				next_value = *it;
				break;
			}
			else if (it->first < min_t) {
				prev_value = *it;
			}
			else {
				values_for_mean.push_back(*it);
			}
		}
3ab32f47   Benjamin Renard   Imf default value... 
108
 
		float value = NAN;
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
 
		if (!values_for_mean.empty()) {
			//Compute mean
			float sum = 0;
			for (std::vector<std::pair<double,float> >::iterator it = values_for_mean.begin(); it != values_for_mean.end(); ++it) {
				sum += it->second;
			}
			value = sum / (float)values_for_mean.size();
		}
		else {
			if (!isNAN(prev_value.first) && !isNAN(next_value.first) && (next_value.first - prev_value.first <= MAX_GAP_SIZE)) {
				//Compute interpolated value
				value = prev_value.second + (time - prev_value.first) / (next_value.first - prev_value.first) * (next_value.second - prev_value.second);
			}
		}

		return value;
	}

	void getImfData(double time, float& b_x, float& b_y, float& b_z) {
3ab32f47   Benjamin Renard   Imf default value... 
128
129
130
 
		b_x = getValue(_b_x_gse,time);
		b_y = getValue(_b_y_gse,time);
		b_z = getValue(_b_z_gse,time);
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
131
132
133
 
	}

	void getPswData(double time, float& p_sw) {
3ab32f47   Benjamin Renard   Imf default value... 
134
 
		p_sw = getValue(_psw,time);
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
135
136
137
 
	}

	void getDstData(double time, float& dst) {
3ab32f47   Benjamin Renard   Imf default value... 
138
 
		dst = getValue(_dst,time);
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
 
	}

private:
	ParamDataSpec<std::vector<float> >& _paramImfInput;

	ParamDataSpec<float>& _paramPswInput;

	ParamDataSpec<float>& _paramDstInput;

	std::vector<std::pair<double,float> > _b_x_gse;

	std::vector<std::pair<double,float> > _b_y_gse;

	std::vector<std::pair<double,float> > _b_z_gse;

	std::vector<std::pair<double,float> > _psw;

	std::vector<std::pair<double,float> > _dst;
};

template <typename ElemType> 
class Tsyganenko96 : public Tsyganenko96Base {
public:
	/**
	 * @brief Constructor.
	 * @details Create the ParamData type of the input ParamData.
	 */
	Tsyganenko96(Process& pProcess, ParamDataSpec<std::vector<ElemType> >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, ParamData& paramDstInput)
	: Tsyganenko96Base(pProcess, paramImfInput, paramPswInput, paramDstInput),
	  _paramInput(paramInput),
	  _paramOutput(new ParamDataSpec<std::vector<ElemType> >) {
	  _paramDataOutput=_paramOutput;
	}

	virtual ~Tsyganenko96() {
	}

	/**
	 * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
	 */
  
	void write(ParamDataIndexInfo &pParamDataIndexInfo) {
		for (unsigned int _index = pParamDataIndexInfo._startIndex ;
			_index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
			++_index)
		{
			double crtTime = _paramInput.getTime(_index);

			float b_x_gse, b_y_gse, b_z_gse;
			getImfData(crtTime, b_x_gse, b_y_gse, b_z_gse);
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
190
191
 
			float p_sw;
			getPswData(crtTime, p_sw);
3ab32f47   Benjamin Renard   Imf default value... 
192
193
194
 
			if (isNAN(p_sw)) {
				p_sw = DEFAULT_PSW;
			}
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
195
196
197
 

			float dst;
			getDstData(crtTime, dst);
3ab32f47   Benjamin Renard   Imf default value... 
198
199
200
 
			if (isNAN(dst)) {
				dst = DEFAULT_DST;
			}
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
 

			std::vector<ElemType> inputElt = _paramInput.get(_index);

			time_t timestamp = crtTime;
			struct tm *tmp;
			tmp = gmtime(&timestamp);

			std::vector<ElemType> ouputElt;
			ouputElt.resize(3);
			ouputElt << NotANumber();

			//Init geopack with GSM frame
			geopack::GeopackWrapper::initInGSM(1900 + tmp->tm_year, 1 + tmp->tm_yday, tmp->tm_hour, tmp->tm_min, tmp->tm_sec);

			//Compute position in GSM frame
			int transform_flag = -1;

			float sat_pos_X_GSE = inputElt[0];
			float sat_pos_Y_GSE = inputElt[1];
			float sat_pos_Z_GSE = inputElt[2];
			float sat_pos_X_GSM = 0.;
			float sat_pos_Y_GSM = 0.;
			float sat_pos_Z_GSM = 0.;

			geogsw_08_(&sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM,
				&sat_pos_X_GSE, &sat_pos_Y_GSE, &sat_pos_Z_GSE, &transform_flag);

			//Check if in magnetopause
			if (geopack::GeopackWrapper::isInMagnetopause(p_sw, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM)) {
				//Compute Imf B field in GSM frame
				float b_x_gsm, b_y_gsm, b_z_gsm;
3ab32f47   Benjamin Renard   Imf default value... 
232
233
234
235
236
237
238
239
240
 
				if (!isNAN(b_x_gse) && !isNAN(b_y_gse) && !isNAN(b_z_gse)) {
					geogsw_08_(&b_x_gsm, &b_y_gsm, &b_z_gsm,
						&b_x_gse, &b_y_gse, &b_z_gse, &transform_flag);
				}
				else {
					b_x_gsm = DEFAULT_IMF_GSM_X;
					b_y_gsm = DEFAULT_IMF_GSM_Y;
					b_z_gsm = DEFAULT_IMF_GSM_Z;
				}
8819b2bc   Benjamin Renard   Add plugin Tsygan... 
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
 

				//Compute magnetic field in GSM
				float B_X_GSM_RES, B_Y_GSM_RES, B_Z_GSM_RES;
				geopack::GeopackWrapper::computeGeomagneticFieldInGSM(inputElt[0], inputElt[1], inputElt[2], p_sw, dst, b_y_gsm, b_z_gsm, 
					B_X_GSM_RES, B_Y_GSM_RES, B_Z_GSM_RES);
				ouputElt[0] = B_X_GSM_RES;
				ouputElt[1] = B_Y_GSM_RES;
				ouputElt[2] = B_Z_GSM_RES;
			}

			_paramOutput->pushTime(crtTime);
			_paramOutput->getDataList().push_back(ouputElt);
		}
	}

private:
	ParamDataSpec<std::vector<ElemType> >& _paramInput;

	ParamDataSpec<std::vector<ElemType> >* _paramOutput;
};  

} /* namespace Tsyganenko96 */
} /* namespace Parameters */
} /* namespace AMDA */
#endif /* TSYGANENKO96_HH_ */