/*
 * 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"
#include "shue_compute.hh"

#include <iterator>
namespace AMDA {
    namespace Parameters {
        namespace Shue {

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

            class ShueBase : public Operation {
            public:

                ShueBase(Process& pProcess, ParamData& paramImfInput, ParamData& paramPswInput)
                : Operation(pProcess),
                _paramImfInput(dynamic_cast<ParamDataSpec<std::vector<float> >&> (paramImfInput)),
                _paramPswInput(dynamic_cast<ParamDataSpec<float>&> (paramPswInput)) {
                }

                virtual ~ShueBase() {
                }

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

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

            private:

                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;

            };

            template <typename ElemType, class TOutputParamData>
            class Shue : public ShueBase {
            public:

                Shue(Process& pProcess, ParamDataSpec<std::vector<ElemType> >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98) :
                ShueBase(pProcess, paramImfInput, paramPswInput),
                _paramInput(paramInput),
                _paramOutput(new TOutputParamData), _inGSE(inGSE), _model98(model98) {
                    _paramDataOutput = _paramOutput;
                }

                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;
                        /**
                         *  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;
                         *   */
                        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);
                        }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);
                        }                        
                        if (computed) {
                            dst = (float) dst*(-pos_id);
                            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;
                            }
                        }
                        _paramOutput->pushTime(crtTime);
                        pushData(ouputElt, dst, pos_id);
                    }
                };

                virtual void pushData(std::vector<ElemType> ouputElt, ElemType dst, int pos_id) = 0;
            protected:
                ParamDataSpec<std::vector<ElemType> >& _paramInput;

                TOutputParamData* _paramOutput;

                bool _inGSE;
                
                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);
                }
            };
            
            
            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);
                }
            };

        } // end Shue
    } // end Parameters
} // end AMDA


#endif /* SHUE_HH */