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

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

				//Compute magnetic field
				float B_X_RES, B_Y_RES, B_Z_RES;
				bool computed = false;
				if (_inGSE) {
					// In GSE
					computed = geopack::GeopackWrapper::computeGeomagneticFieldInGSE(inputElt[0], inputElt[1], inputElt[2], 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(inputElt[0], inputElt[1], inputElt[2], 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;
};  

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