/*
 * 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)

#define DEFAULT_IMF_GSM_X 0.
#define DEFAULT_IMF_GSM_Y 2.
#define DEFAULT_IMF_GSM_Z -3.

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

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

                void getDstData(double time, float& dst) {
                    dst = getValue(_dst, time);
                }

            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, bool inGSE)
                : Tsyganenko96Base(pProcess, paramImfInput, paramPswInput, paramDstInput),
                _paramInput(paramInput),
                _paramOutput(new ParamDataSpec<std::vector<ElemType> >), _inGSE(inGSE) {
                    _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);

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

                        float dst;
                        getDstData(crtTime, dst);
                        if (isNAN(dst)) {
                            dst = DEFAULT_DST;
                        }

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

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

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

                            //Compute magnetic field
                            float B_X_RES, B_Y_RES, B_Z_RES;
                            bool computed = false;
                            if (_inGSE) {
                                // In GSE
                                computed = geopack::GeopackWrapper::computeGeomagneticFieldInGSE(sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM, p_sw, dst, b_y_gsm, b_z_gsm,
                                        B_X_RES, B_Y_RES, B_Z_RES);
                            } else {
                                // In GSM
                                computed = geopack::GeopackWrapper::computeGeomagneticFieldInGSM(sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM, p_sw, dst, b_y_gsm, b_z_gsm,
                                        B_X_RES, B_Y_RES, B_Z_RES);
                            }
                            if (computed) {
                                ouputElt[0] = B_X_RES;
                                ouputElt[1] = B_Y_RES;
                                ouputElt[2] = B_Z_RES;
                            }
                        }

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

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

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

                bool _inGSE;
            };

            template<typename ElemType, class TOutputParamData>

            class DipoleBase : public Operation {
            public:

                DipoleBase(Process& pProcess, ParamDataSpec<ElemType>& paramInput) : Operation(pProcess),

                _paramInput(paramInput),

                _paramOutput(new TOutputParamData) {

                    _paramDataOutput = _paramOutput;
                }

                virtual ~DipoleBase() {
                }

                void write(ParamDataIndexInfo &pParamDataIndexInfo) {
                    for (unsigned int _index = pParamDataIndexInfo._startIndex;
                            _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
                            ++_index) {
                        double crtTime = _paramInput.getTime(_index);
                        time_t timestamp = crtTime;
                        struct tm *tmp;
                        tmp = gmtime(&timestamp);
                        float tiltAngle;
                        std::vector<float> attitudeGSE;
                        attitudeGSE.resize(3);
                        std::vector<float> attitudeGSM;
                        attitudeGSM.resize(3);
                        geopack::GeopackWrapper::computeTiltAngleAttitude(1900 + tmp->tm_year, 1 + tmp->tm_yday, tmp->tm_hour, tmp->tm_min, tmp->tm_sec,
                                tiltAngle,
                                attitudeGSE,
                                attitudeGSM);
                        ElemType tiltAngle_res = (ElemType) tiltAngle;
                        std::vector<ElemType> attitudeGSE_res;
                        attitudeGSE_res.resize(3);
                        attitudeGSE_res[0] = (ElemType) attitudeGSE[0];
                        attitudeGSE_res[1] = (ElemType) attitudeGSE[1];
                        attitudeGSE_res[2] = (ElemType) attitudeGSE[2];

                        std::vector<ElemType> attitudeGSM_res;
                        attitudeGSM_res.resize(3);
                        attitudeGSM_res[0] = (ElemType) attitudeGSM[0];
                        attitudeGSM_res[1] = (ElemType) attitudeGSM[1];
                        attitudeGSM_res[2] = (ElemType) attitudeGSM[2];

                        attitudeGSM_res.resize(3);
                        _paramOutput->pushTime(crtTime);
                        pushData(tiltAngle_res, attitudeGSE_res, attitudeGSM_res);
                    }
                };
                virtual void pushData(ElemType tilt, std::vector<ElemType> attitudeGSE, std::vector<ElemType> attitudeGSM) = 0;

            protected:
                /**
                 * @brief Input paramter data.
                 */

                ParamDataSpec<ElemType>& _paramInput;

                TOutputParamData* _paramOutput;

            };

            template <typename ElemType>

            class DipoleTiltAngle : public DipoleBase<ElemType, ParamDataSpec<ElemType>> {

            public:

                /**
                 * @brief Constructor.
                 * @details Create the ParamData type of the input ParamData.
                 */

                DipoleTiltAngle(Process& pProcess, ParamDataSpec<ElemType >& paramInput)

                : DipoleBase<ElemType , ParamDataSpec<ElemType> >(pProcess, paramInput) {

                }

                virtual ~DipoleTiltAngle() {

                }

                void pushData(ElemType tilt, std::vector<ElemType> /*attitudeGSE*/, std::vector<ElemType> /*attitudeGSM*/) {
                    DipoleBase<ElemType, ParamDataSpec<ElemType>>::_paramOutput->getDataList().push_back(tilt * 180.0 / (atan(1)*4));
                }
            };
            
           template <typename ElemType>
            class DipoleAtitudeGSM : public DipoleBase<ElemType, ParamDataSpec<std::vector<ElemType>>> {
            public:

                /**
                 * @brief Constructor.
                 * @details Create the ParamData type of the input ParamData.
                 */
                DipoleAtitudeGSM(Process& pProcess, ParamDataSpec<ElemType >& paramInput)
                : DipoleBase<ElemType, ParamDataSpec<std::vector<ElemType>>>(pProcess, paramInput) {
                }

                virtual ~DipoleAtitudeGSM() {
                }

                void pushData(ElemType /*tilt*/, std::vector<ElemType> /*attitudeGSE*/, std::vector<ElemType> attitudeGSM) {
                    DipoleBase<ElemType, ParamDataSpec<std::vector<ElemType>>>::_paramOutput->getDataList().push_back(attitudeGSM);
                }
            };
            
                       template <typename ElemType>
            class DipoleAtitudeGSE : public DipoleBase<ElemType, ParamDataSpec<std::vector<ElemType>>> {
            public:

                /**
                 * @brief Constructor.
                 * @details Create the ParamData type of the input ParamData.
                 */
                DipoleAtitudeGSE(Process& pProcess, ParamDataSpec<ElemType >& paramInput)
                : DipoleBase<ElemType, ParamDataSpec<std::vector<ElemType>>>(pProcess, paramInput) {
                }

                virtual ~DipoleAtitudeGSE() {
                }

                void pushData(ElemType /*tilt*/, std::vector<ElemType> attitudeGSE, std::vector<ElemType> /*attitudeGSM*/) {
                    DipoleBase<ElemType , ParamDataSpec<std::vector<ElemType>>>::_paramOutput->getDataList().push_back(attitudeGSE);
                }

            };

        } /* namespace Tsyganenko96 */
    } /* namespace Parameters */

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