CDPP / AMDA_Kernel

Blame view

src/ExternLib/Shue/Shue.hh 13 KB
Edit Raw Normal View History
45fb3489   Hacene SI HADJ MOHAND   plugin 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/* 
 * File:   Shue.hh
 * Author: hacene
 *
 * Created on December 31, 2020, 10:14 AM
 */

#ifndef SHUE_HH
#define SHUE_HH

#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
21
 
#include "shue_compute.hh"
45fb3489   Hacene SI HADJ MOHAND   plugin 
22
23
24
25
 

#include <iterator>
namespace AMDA {
    namespace Parameters {
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
26
27
 
        namespace Shue {
45fb3489   Hacene SI HADJ MOHAND   plugin 
28
29
30
31
32
33
 
#define AVERAGE_TIME 1200 // (seconds)
#define MAX_GAP_SIZE 3600 // (seconds)
#define DEFAULT_IMF_GSM_X 0.
#define DEFAULT_IMF_GSM_Y 2.
#define DEFAULT_IMF_GSM_Z -3.
#define DEFAULT_PSW 3.
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
34
35
 

            class ShueBase : public Operation {
45fb3489   Hacene SI HADJ MOHAND   plugin 
36
 
            public:
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
37
45fb3489   Hacene SI HADJ MOHAND   plugin 
38
39
40
 
                ShueBase(Process& pProcess, ParamData& paramImfInput, ParamData& paramPswInput)
                : Operation(pProcess),
                _paramImfInput(dynamic_cast<ParamDataSpec<std::vector<float> >&> (paramImfInput)),
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
41
 
                _paramPswInput(dynamic_cast<ParamDataSpec<float>&> (paramPswInput)) {
45fb3489   Hacene SI HADJ MOHAND   plugin 
42
 
                }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
43
44
 

                virtual ~ShueBase() {
45fb3489   Hacene SI HADJ MOHAND   plugin 
45
 
                }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
46
45fb3489   Hacene SI HADJ MOHAND   plugin 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
 
                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));
                    }
                }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
70
45fb3489   Hacene SI HADJ MOHAND   plugin 
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
 
                float getValue(std::vector<std::pair<double, float> >& input, double time) {
                    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);
                        }
                    }

                    float value = NAN;
                    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) {
                    b_x = getValue(_b_x_gse, time);
                    b_y = getValue(_b_y_gse, time);
                    b_z = getValue(_b_z_gse, time);
                }

                void getPswData(double time, float& p_sw) {
                    p_sw = getValue(_psw, time);
                }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
117
118
119
 

            private:
45fb3489   Hacene SI HADJ MOHAND   plugin 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
 
                ParamDataSpec<std::vector<float> >& _paramImfInput;

                ParamDataSpec<float>& _paramPswInput;


                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;

            };
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
134
135
136
 

            template <typename ElemType, class TOutputParamData>
            class Shue : public ShueBase {
45fb3489   Hacene SI HADJ MOHAND   plugin 
137
 
            public:
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
138
139
 

                Shue(Process& pProcess, ParamDataSpec<std::vector<ElemType> >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98) :
45fb3489   Hacene SI HADJ MOHAND   plugin 
140
 
                ShueBase(pProcess, paramImfInput, paramPswInput),
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
141
142
 
                _paramInput(paramInput),
                _paramOutput(new TOutputParamData), _inGSE(inGSE), _model98(model98) {
45fb3489   Hacene SI HADJ MOHAND   plugin 
143
144
 
                    _paramDataOutput = _paramOutput;
                }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
 

                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);

                        float p_sw;
                        getPswData(crtTime, p_sw);
                        if (isNAN(p_sw)) {
                            p_sw = DEFAULT_PSW;
                        }


                        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
                        ShueCompute::initInGSM(1900 + tmp->tm_year, 1 + tmp->tm_yday, tmp->tm_hour, tmp->tm_min, tmp->tm_sec);

                        //Compute shue 98 
                        // inout intialisation

                        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.;

                        gswgse_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);

                        //Compute Imf B field in GSM frame
                        float b_x_gsm, b_y_gsm, b_z_gsm;
                        if (!isNAN(b_x_gse) && !isNAN(b_y_gse) && !isNAN(b_z_gse)) {
                            gswgse_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;
                        }

                        float X_GSW, Y_GSW, Z_GSW;
                        int  flag_transToGSE = 1;
                        float X_GSE = 0;
                        float Y_GSE = 0;
                        float Z_GSE = 0;
                        
                        float dst;
                        int pos_id;
                        bool computed = false;
2cae24f0   Hacene SI HADJ MOHAND   shue 97 ok 
212
213
214
215
216
217
 
                        /**
                         *  for testing to reproduce JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 103, NO. A8, PAGES 17,691-17,700, AUGUST 1, 1998
                         *  fig 4         
                        p_sw=14.0; 
                        b_z_gsm = 9.3;
                         *   */
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
218
219
220
 
                        if(_model98){
                        computed = ShueCompute::shue98(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM,
                                X_GSW, Y_GSW, Z_GSW, dst, pos_id);
2cae24f0   Hacene SI HADJ MOHAND   shue 97 ok 
221
222
223
224
 
                        }else{
                                 computed = ShueCompute::shue97(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM,
                                X_GSW, Y_GSW, Z_GSW, dst, pos_id);
                        }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
225
 
                        if (computed) {
9582722d   Hacene SI HADJ MOHAND   adding sign to th... 
226
 
                            dst = (float) dst*(-pos_id);
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
227
228
229
230
231
232
233
234
235
236
237
238
239
 
                            if(_inGSE){
                             // transform results from GSM to GSE
                            gswgse_08_(&X_GSW, &Y_GSW, &Z_GSW, 
                                                  &X_GSE, &Y_GSE,   &Z_GSE, &flag_transToGSE);
                            ouputElt[0] = X_GSE;
                            ouputElt[1] = Y_GSE;
                            ouputElt[2] = Z_GSE;    
                            }else{
                            ouputElt[0] = X_GSW;
                            ouputElt[1] = Y_GSW;
                            ouputElt[2] = Z_GSW;
                            }
                        }
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
240
241
242
243
244
245
246
 
                        _paramOutput->pushTime(crtTime);
                        pushData(ouputElt, dst, pos_id);
                    }
                };

                virtual void pushData(std::vector<ElemType> ouputElt, ElemType dst, int pos_id) = 0;
            protected:
45fb3489   Hacene SI HADJ MOHAND   plugin 
247
248
 
                ParamDataSpec<std::vector<ElemType> >& _paramInput;
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
249
 
                TOutputParamData* _paramOutput;
45fb3489   Hacene SI HADJ MOHAND   plugin 
250
251
 

                bool _inGSE;
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
 
                
                bool _model98;

            };

            template <typename ElemType>
            class ShuePos : public Shue<ElemType, ParamDataSpec<std::vector<ElemType>>>
            {
            public:

                ShuePos(Process& pProcess, ParamDataSpec<std::vector<ElemType> >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): 
                Shue<ElemType, ParamDataSpec<std::vector < ElemType>>>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){
                }
                virtual ~ShuePos(){
                    
                }
              
                void pushData(std::vector<ElemType> ouputElt, ElemType /*dst*/, int /*pos_id*/) {
                    Shue<ElemType, ParamDataSpec<std::vector < ElemType>>>::_paramOutput->getDataList().push_back(ouputElt);
                }
45fb3489   Hacene SI HADJ MOHAND   plugin 
272
 
            };
45fb3489   Hacene SI HADJ MOHAND   plugin 
273
274
5d3b0cd0   Hacene SI HADJ MOHAND   shue98 ok 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
 
            template <typename ElemType>
            class ShueDst : public Shue<ElemType, ParamDataSpec<ElemType>>
            {
            public:

                ShueDst(Process& pProcess, ParamDataSpec<std::vector<ElemType> >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): 
                Shue<ElemType, ParamDataSpec<ElemType>>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){
                }
                virtual ~ShueDst(){
                    
                }
              
                void pushData(std::vector<ElemType> /*ouputElt*/, ElemType dst, int /*pos_id*/) {
                    Shue<ElemType, ParamDataSpec< ElemType>>::_paramOutput->getDataList().push_back(dst);
                }
            };
            
           template <typename ElemType>
            class ShuePos_flag : public Shue<ElemType, ParamDataSpec<int>>
            {
            public:

               ShuePos_flag(Process& pProcess, ParamDataSpec<std::vector<ElemType> >&  paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): 
                Shue<ElemType, ParamDataSpec<int>>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){
                }
                virtual ~ShuePos_flag(){
                    
                }
              
                void pushData(std::vector<ElemType> /*ouputElt*/, ElemType /*dst*/, int pos_id) {
                    Shue<ElemType, ParamDataSpec<int>>::_paramOutput->getDataList().push_back(pos_id);
                }
            };
45fb3489   Hacene SI HADJ MOHAND   plugin 
309
310
311
312
313
314
 
        } // end Shue
    } // end Parameters
} // end AMDA


#endif /* SHUE_HH */