/*
 * Morschhauser.hh
 *
 *  Created on: May 30, 2016
 *      Author: AKKA IS
 */

#ifndef MORSCHHAUSER_HH_
#define MORSCHHAUSER_HH_

#include "morschhauser_constants.hh"

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

namespace AMDA {
namespace Parameters {

/**
 * @class MorschhauserCommon
 * @brief It is responsible to compute Morshhauser Mars Model Magnetic Field along an orbit. Abstract class
 * @details This class implement the interface Operation. Input vector given in MSO coordinates system.
 */
template<typename DataType, class TOutputParamData>
class MorschhauserCommon : public Operation {
public:
	/**
	 * @brief Constructor.
	 * @details Create the ParamData type of the input ParamData.
	 */
	MorschhauserCommon(Process& pProcess, ParamDataSpec<std::vector<DataType> >&paramInput): Operation(pProcess),
	  _paramInput(paramInput),
	  _paramOutput(new TOutputParamData())
	{
		_paramDataOutput=_paramOutput;

		init_coef_V_int();
		init_coef_V_ext_night();
		init_coef_V_ext_day();
	}

	virtual ~MorschhauserCommon() {
	}

	/**
	 * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
	 */
	void write(ParamDataIndexInfo &pParamDataIndexInfo) {
		for (unsigned int _index = pParamDataIndexInfo._startIndex ;
				_index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
				++_index)
		{
			std::vector<DataType> lVal = _paramInput.getDataList()[_index];
			short sign = (lVal[0] > 0) ? 1 : 0;

			//BRE - ToDo Corrdinates transformation before magnetic field calculation
			std::vector<DataType> outputB;
			outputB.resize(3);

			DataType outputBm;
			field_(lVal, sign, outputB, outputBm);

			_paramOutput->pushTime(_paramInput.getTime(_index));
			pushResult(outputB, outputBm);
		}
	}

protected:

	virtual void pushResult(std::vector<DataType> field, DataType magnitude) = 0;

	/**
	 * @brief Input paramter data.
	 */
	ParamDataSpec<std::vector<DataType> >& _paramInput;

	/**
	 * @brief Output parameter data.
	 */
	TOutputParamData *_paramOutput;

private:
	long double _coef_V_int[111][221], _coef_V_ext_day[6][11], _coef_V_ext_night[11][21];

	int field_(std::vector<DataType> inputElt, short sign, std::vector<DataType>& outputB, DataType& outputBm) {
		long double r,dtheta,dphi,x,y,z;

		long double cos_theta, sin_theta, cos_phi, sin_phi;

		long double P[111][111]={ {0.0L} , {0.0L} };
		long double Br, Btheta, Bphi;

		x =  (long double)inputElt[0];
		y =  (long double)inputElt[1];
		z =  (long double)inputElt[2];

		r = sqrtl( x * x + y * y + z * z );

		dtheta = acosl( z / r );
		dphi   = atanl( y / x );

		if ( x < 0.0 )
		{
			dphi += M_PIl;
		}

		if ( r < Rlim )
		{
			cos_theta = cosl( dtheta );
			sin_theta = sinl( dtheta );
			cos_phi = cosl( dphi );
			sin_phi = sinl( dphi );

			plgndr( cos_theta, P );
			V_int( Rma * r , dtheta , dphi, sign, P, Br, Btheta, Bphi);
			outputBm   = (double)sqrtl( Br * Br + Btheta * Btheta + Bphi * Bphi );
			outputB.push_back((double)(Br * cos_phi * sin_theta - Bphi * sin_phi  + Btheta * cos_phi * cos_theta));
			outputB.push_back( (double)(Br * sin_phi * sin_theta + Bphi * cos_phi));
			outputB.push_back((double)(Br * cos_theta - Btheta * sin_theta));
		}
		else
		{
			Br = ( Btheta = ( Bphi = 0.0L / 0.0L ) );
			outputBm << NotANumber();
			outputB << NotANumber();
		}
		return 0;
	}

	void init_coef_V_int(void) {
		for (int i = 0; i < V_INT_SIZE; ++i)
		{
			_coef_V_int[(int)V_INT[i][0]][(int)(V_INT[i][0]+V_INT[i][1])] = V_INT[i][2];
		}
	}

	void init_coef_V_ext_night(void) {
		for (int i = 0; i < V_EXT_NIGHT_SIZE; ++i)
		{
			_coef_V_ext_night[(int)V_EXT_NIGHT[i][0]][(int)(V_EXT_NIGHT[i][0]+V_EXT_NIGHT[i][1])] = V_EXT_NIGHT[i][2];
		}
	}

	void init_coef_V_ext_day(void) {
		for (int i = 0; i < V_EXT_DAY_SIZE; ++i)
		{
			_coef_V_ext_day[(int)V_EXT_DAY[i][0]][(int)(V_EXT_DAY[i][0]+V_EXT_DAY[i][1])] = V_EXT_DAY[i][2];
		}
	}

	int  V_int(long double r , long double theta , long double phi, short sign, long double P[111][111],
			long double& Br, long double& Btheta, long double& Bphi) {
		long double V1,V2,V3,V4,V5,V6,x,y,z,Ypm,Ymm, puis,Ypmav,Ymmav,puis2;
		int l,m;
		long double tan_theta, sin_theta, cos_theta;

		tan_theta = tanl( theta );
		sin_theta = sinl( theta );
		cos_theta = cosl( theta );

		Br     = 0.0L;
		Bphi   = 0.0L;
		Btheta = 0.0L;
		x = ( y = ( z = 0.0L ) );
		Ypm = ( Ypmav = 0.0L );
		Ymm = ( Ymmav = 0.0L );

		for (l = 1 ; l <= 110 ; l++ )
		{
			V1 = 0.0L;
			V2 = 0.0L;
			V3 = 0.0L;
			V4 = 0.0L;
			V5 = 0.0L;
			V6 = 0.0L;
			puis = 0.0L;
			Ypmav = 0.0L;
			Ymmav = 0.0L;

			for ( m = l ; m > 0 ; m-- )
			{
				Ypm = Ylm( l ,  m , theta , phi, P );
				Ymm = Ylm( l , -m , theta , phi, P );

				V1 += _coef_V_int[l][l+m] * Ypm;
				V1 += _coef_V_int[l][l-m] * Ymm;

				V2 -= m * _coef_V_int[l][l+m] * Ymm;
				V2 += m * _coef_V_int[l][l-m] * Ypm;

				if( sign == 0 && l <= 10 ) // night side MSO
				{
					V4 += _coef_V_ext_night[l][l+m] * Ypm;
					V4 += _coef_V_ext_night[l][l-m] * Ymm;
					V5 -= m * _coef_V_ext_night[l][l+m] * Ymm;
					V5 += m * _coef_V_ext_night[l][l-m] * Ypm;
					V6 += _coef_V_ext_night[l][l+m] * ( m * Ypm / tan_theta - sqrtl( l - m + l * l - m * m ) * cosl( m * phi ) / cosl( ( m + 1.0l ) * phi ) * Ypmav );
					V6 += _coef_V_ext_night[l][l-m] * ( m * Ymm / tan_theta - sqrtl( l - m + l * l - m * m ) * sinl( m * phi ) / sinl( ( m + 1.0l ) * phi ) * Ymmav );
				}

				if( sign > 0 && l <= 5 )  // day side MSO
				{
					V4 += _coef_V_ext_day[l][l+m] * Ypm;
					V4 += _coef_V_ext_day[l][l-m] * Ymm;
					V5 -= m * _coef_V_ext_day[l][l+m] * Ymm;
					V5 += m * _coef_V_ext_day[l][l-m] * Ypm;
					V6 += _coef_V_ext_day[l][l+m] * ( m * Ypm / tan_theta - sqrtl( l - m + l * l - m * m ) * cosl( m * phi ) / cosl( ( m + 1.0l ) * phi ) * Ypmav );
					V6 += _coef_V_ext_day[l][l-m] * ( m * Ymm / tan_theta - sqrtl( l - m + l * l - m * m ) * sinl( m * phi ) / sinl( ( m + 1.0l ) * phi ) * Ymmav );
				}

				if ( phi != 0.0L )
				{
					V3 += _coef_V_int[l][l+m] * ( m * Ypm / tan_theta - sqrtl( l - m + l * l - m * m ) * cosl( m * phi ) / cosl( ( m + 1.0l ) * phi ) * Ypmav );
					V3 += _coef_V_int[l][l-m] * ( m * Ymm / tan_theta - sqrtl( l - m + l * l - m * m ) * sinl( m * phi ) / sinl( ( m + 1.0l ) * phi ) * Ymmav );
				}
				else
				{
					V3 += _coef_V_int[l][l+m] * ( m * Ypm / tan_theta - sqrtl( l - m + l * l - m * m ) * Ypmav );
				}

				Ypmav = Ypm;
				Ymmav = Ymm;
			}

			V1 +=   _coef_V_int[l][l] * Ylm( l , 0 , theta , phi, P );

			if( sign == 0 && l <= 10 ) // night side MSO
			{
				V4 +=   _coef_V_ext_night[l][l] * Ylm( l , 0 , theta , phi, P );
			}

			if( sign > 0 && l <= 5 ) // day side MSO
			{
				V4 +=   _coef_V_ext_day[l][l] * Ylm( l , 0 , theta , phi, P );
			}

			V3 +=   _coef_V_int[l][l] * l / sin_theta * ( cos_theta * P[l][0]-P[l-1][0] );

			if( sign == 0 && l <= 10 ) // night side MSO
			{
				V6 +=   _coef_V_ext_night[l][l] * l / sin_theta * ( cos_theta * P[l][0]-P[l-1][0] );
			}

			if( sign > 0 && l <= 5 ) // day side MSO
			{
				V6 +=   _coef_V_ext_day[l][l] * l / sin_theta * ( cos_theta * P[l][0]-P[l-1][0] );
			}


			puis = powl( Rma / r , l + 2.0L );
			puis2= powl( r / Rma , l - 1.0L );

			z -= - ( l + 1.0L ) * V1 * puis + l * V4 * puis2; //minus sign disappear because magnetic field is the ooposite of the gradient
			x -=                  ( V2 * puis + V5 * puis2 ) / sin_theta;
			y -=                  ( V3 * puis + V6 * puis2 );

		}

		Br = z;
		Bphi = x;
		Btheta= y;

		return 0;
	}

	void plgndr(double x, long double P[111][111]) {
		int m,l;
		double fact,coef;

		fact = 1.0;
		coef = sqrt( 1.0 - x * x );
		P[0][0] = 1.0;


		for ( m = 1 ; m <= 110 ; m++ )
		{

			P[m][m] = - P[m-1][m-1] * fact * coef;
			fact += 2.0;
		}

		for ( m = 0 ; m <= 109 ; m++ )
		{
			P[m+1][m] = x * ( 2.0 * m + 1.0 ) * P[m][m];
		}

		for ( m = 0 ; m <= 108 ; m++ )
		{
	        for ( l = m + 2 ; l <= 110 ; l++ )
	        {
	        	P[l][m] = ( x * ( 2.0 * l - 1.0 ) * P[l-1][m] - ( l + m - 1.0 ) * P[l-2][m] ) / ( l - m );
	        }

		}
	}

	long double Ylm(int l , int m , long double /*theta*/ , long double phi, long double P[111][111]) {
		long double Plm;
		int k;
		long double p, sqrt2 = sqrtl( 2.0L );

		Plm = P[l][abs(m)];

		p = 1;

		if ( m != 0 )
		{
	        for ( k = (l - abs(m) + 1 ) ; k <= ( l + abs(m) ) ; k++)
	        {
	        	p *= sqrtl( (long double) k );
	        }
	        Plm /= p;
	        Plm *= sqrt2;
		}

		if ( m < 0 )
		{
			Plm *= sinl( fabsl( m ) * phi );
		}
		else if ( m > 0 )
		{
	        Plm *= cosl(  fabsl( m ) * phi );
		}

		Plm *= powl( - 1.0L , fabsl( m ) );

		if ( abs(m) > l )
		{
			Plm = 0.0L;
		}

		return Plm;
	}
};

/**
 * @class MorschhauserBField
 * @brief It is responsible to compute Morshhauser Mars Model Magnetic Field along an orbit.
 * @details Input vector given in MSO coordinates system.
 */
template<typename DataType>
class MorschhauserBField : public MorschhauserCommon<DataType, ParamDataSpec<std::vector<DataType> >> {
public:
	/**
	 * @brief Constructor.
	 * @details Create the ParamData type of the input ParamData.
	 */
	MorschhauserBField(Process& pProcess, ParamDataSpec<std::vector<DataType> >& paramInput) : MorschhauserCommon<DataType, ParamDataSpec<std::vector<DataType> >>(pProcess, paramInput) {
	}

	virtual ~MorschhauserBField() {

	}

protected:
	virtual void pushResult(std::vector<DataType> field, DataType /*magnitude*/) {
		MorschhauserCommon<DataType, ParamDataSpec<std::vector<DataType> >>::_paramOutput->getDataList().push_back(field);
	}
};

/**
 * @class MorschhauserBMag
 * @brief It is responsible to compute Morshhauser Mars Model Magnetic Field along an orbit.
 * @details Input vector given in MSO coordinates system.
 */
template<typename DataType>
class MorschhauserBMag : public MorschhauserCommon<DataType, ParamDataSpec<DataType>> {
public:
	/**
	 * @brief Constructor.
	 * @details Create the ParamData type of the input ParamData.
	 */
	MorschhauserBMag(Process& pProcess, ParamDataSpec<std::vector<DataType> >& paramInput) : MorschhauserCommon<DataType, ParamDataSpec<DataType>>(pProcess, paramInput) {
	}

	virtual ~MorschhauserBMag() {
	}

protected:
	virtual void pushResult(std::vector<DataType> /*field*/, DataType magnitude) {
		MorschhauserCommon<DataType, ParamDataSpec<DataType>>::_paramOutput->getDataList().push_back(magnitude);
	}
};

} /* namespace Parameters */
} /* namespace AMDA */
#endif /* MORSCHHAUSER_HH_ */