/*
 * 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:   getjupiterMag.hh
 * Author: hacene
 *
 * Created on May 29, 2020, 2:50 PM
 */

#ifndef Jupiter_JRM09_CAN81_HH
#define Jupiter_JRM09_CAN81_HH

#include "Jupiter_JRM09_CAN81_tools.hh"

#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"
#include "SpiceKernelMgr.hh"
#include <vector>

using namespace std;

namespace AMDA {
    namespace Parameters {

        /**
         * @class GetjupiterMagCommon
         * @brief It is responsible to compute Jupiter Magnetic Field along an orbit. Abstract class
         * @details This class implement the interface Operation. Input vector are [R, Z] 
         */

        template<typename DataType, class TOutputParamData>
        class Jupiter_JRM09_CAN81 : public Operation {
        public:

            /**
             * @brief Constructor.
             * @details Create the ParamData type of the input ParamData.
             */
            Jupiter_JRM09_CAN81(Process& pProcess, ParamDataSpec<vector<DataType> >&paramInput, bool can) : Operation(pProcess),
            _paramInput(paramInput),
            _paramOutput(new TOutputParamData()), _can(can) {
                _paramDataOutput = _paramOutput;
            }

            virtual ~Jupiter_JRM09_CAN81() {
            }

            /**
             * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
             */
            void write(ParamDataIndexInfo &pParamDataIndexInfo) {
                AMDA::GETjupiterMAG_TOOLS::jupiter_JRM09_CAN81_tools tools;
                
                for (unsigned int _index = pParamDataIndexInfo._startIndex;
                        _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
                        ++_index) {
                    vector<DataType> in = _paramInput.getDataList()[_index];
                    if (in[0] == 0)
                        in[0] = 1e-31;
                    vector<DataType> lVal;
                    lVal.resize(3);
                    lVal[0] = std::sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]);
                    if (lVal[0] == 0)
                        lVal[0] = 1e-31;
                    // 
                    lVal[1] = std::acos(in[2] / lVal[0])*180 / PI;
                    //
                    
                    lVal[2] = std::atan2(in[1] , in[0])*180 / PI; 
//                     if (in[0] < 0) 
//                         lVal[2] += 180.0;
                    if (lVal[2] < 0)
                        lVal[2] += 360.0;
                             
                    double crtTime = _paramInput.getTime(_index);

                    vector<DataType> outputBrtp;
                    outputBrtp.resize(3);

                    vector<DataType>outputBxyz;
                    outputBxyz.resize(3);

                    DataType outputBm;

                    // GetJupiterMag00(vector<DataType> inputSphCoord, bool CAN, vector<DataType>& outputBrtp, DataType& outputBm, vector<DataType>& outputBxyz) 
                    GetJupiterMag00(lVal, _can, outputBrtp, outputBm, outputBxyz);
                    _paramOutput->pushTime(crtTime);
                    pushResult(outputBrtp, outputBm, outputBxyz);
                }
            }
        protected:

            virtual void pushResult(vector<DataType> fieldSph, DataType magnitude_, vector<DataType> fieldCoor) = 0;

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

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

            /**
             * use can model or not
             */
            bool _can;

        private:

            int can81(vector<DataType>& inputCoor, vector<DataType>& outpField) {
                /**
                 * @brief Jupiter's external magnetic field model:
                 * @details ;   The Connerney+1981_JGR model: Modeling the Jovian Current Sheet and Inner Magnetophere.
                 *The CAN model should be only used near the Jupiter (<~30 RJ).
                 *Input: rho and z in cylindrical coords, z is along the magnetic axis.
                 *Output: Brho and Bz (nT) also in cylindrical coords. Note that the model is axis-symmetric.
                 */

                //   [constants]
                double D = 2.5; // Rj
                double a = 5.0; // Rj
                double b = 50.0; // 70.0         ; Rj         ; 50.0 in Connerney1981
                double mu0_I0 = 450.0; //;   175.0*2       ; 225.0*2 in Connerney1981

                double rho = inputCoor[0];
                double z = inputCoor[1];
                double F1, F2, F3, tmp, Bz, Brho, Brho2, Bz2;


                // [Approximate formulas given in Connerney+1981_The magnetic field in Jupiter.]
                // [Cylindrical coordinates: in nT]

                if (rho < a) {
                    F1 = sqrt((z - D)*(z - D) + a * a);
                    F2 = sqrt((z + D)*(z + D) + a * a);
                    F3 = sqrt(z * z + a * a);

                    Brho = 0.5 * rho * (1.0 / F1 - 1.0 / F2);
                    tmp = (z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2);
                    Bz = 2.0 * D / F3 - 0.25 * rho * rho*tmp;
                } else {
                    F1 = sqrt((z - D)*(z - D) + rho * rho);
                    F2 = sqrt((z + D)*(z + D) + rho * rho);
                    F3 = sqrt(z * z + rho * rho);

                    Brho = (F1 - F2 + 2 * D) / rho;
                    if (abs(z) > D && z < 0) Brho = (F1 - F2 - 2 * D) / rho;
                    if (abs(z) < D) Brho = (F1 - F2 + 2 * z) / rho;
                    Brho = Brho - 0.25 * a * a * rho * (1.0 / (F1 * F1 * F1) - 1.0 / (F2 * F2 * F2));
                    tmp = (z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2);
                    Bz = 2.0 * D / F3 - 0.25 * a * a*tmp;
                }
                F1 = sqrt((z - D)*(z - D) + b * b);
                F2 = sqrt((z + D)*(z + D) + b * b);
                F3 = sqrt(z * z + b * b);
                Brho2 = 0.5 * rho * (1 / F1 - 1 / F2);
                Bz2 = 2 * D / F3 - 0.25 * rho * rho * ((z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2));

                Brho -= Brho2;
                Bz -= Bz2;

                Brho *= 0.5 * mu0_I0;
                Bz *= 0.5 * mu0_I0;
                outpField[0] = Brho;
                outpField[1] = Bz;

                return 0;
            };

            /**
             * GetJupiterMag00.pro
              =================================================
             Description:
                   This program is to calculate the Jovian magnetic field of JRM09 (Connerney et al., 2018) and
                    CAN (Connerney et al., 1981) model.

                    All of the relevant subroutines are included in this single file, including the Schmdt coefficients
                    needed in JRM09 model.

                   (1) Jupiter's Magnetic field Model: JRM09 (+CAN), CAN.
   (2) Refer to: Connerney+2018_GRL  JRM09 model
                       Connerney+1981            CAN model
   (3) All the coordinates used here are Right-Hand-Side, because of SPICE convention.
   (4) The positions can be in Arrays, i.e., multi-positions at one time, which can be more efficient.

 Input:
   r (Rj), theta (deg), phi (deg), spherical coordinates in System III RHS (i.e., IAU_Jupiter), better use double!

   optional:
       /CAN              include the CAN model or not (true = yes)
 Output:
   Br (nT), Bt (nT), Bp (nT), Bmag (nT) in System III RHS.
   Bxyz = Bxyz   optional output Bxyz in System III RHS

Example:
   GetJupiterMag00, [20.0d0,20.00000001d0], [90.0d0,90.000001d0], [120.0d0, 120.0d0], $
     Br, Bt, Bp, Bmag1, Bxyz = Bxyz, /can
   GetJupiterMag00, 20.0d0, 90.000001d0, 120.0d0, Br, Bt, Bp, Bmag2, Bxyz = Bxyz, /can

 History:
   (1) writen by Yuxian Wang 2020-01-16 02:26:49
   ywang@irap.omp.eu + yxwang@spaceweather.ac.cn
 =================================================
             */

            int GetJupiterMag00(vector<DataType> inputSphCoord, bool CAN, vector<DataType>& outputBrtp, DataType& outputBm, vector<DataType>& outputBxyz) {
                const AMDA::GETjupiterMAG_TOOLS::COEF _COEF;
                AMDA::GETjupiterMAG_TOOLS::jupiter_JRM09_CAN81_tools tools;
                // [position in: ['Rj, deg, deg'] to ['Rj, rad,rad']]
                long double gsmdtlm, hsmdtlm, phifactor, Pnm, dpnmdtheta, sintheta;
                DataType r, theta, phi, x, radialfactor;
                r = inputSphCoord[0];
                theta = inputSphCoord[1] * PI / 180;
                phi = inputSphCoord[2] * PI / 180;
                std::vector<DataType> inputSphCoord_rad {r, theta, phi};

                // [try to avoid invoke restore one more time.]
                x = cos(theta);
                int maxN = 10; //  JRM09 degrees
                vector<vector<DataType> > Plm(maxN + 1, vector<DataType>(maxN + 1, 0));
                vector<vector < DataType >> dPlm(maxN + 1, vector<DataType>(maxN + 1, 0));

                //int AllLegendre(T1 x, int lma, int mmax, vector<vector<T2>> Plm,
                //vector<vector<T2>> dPlm, bool isSineTheta, int SineTheta, bool DoDerivative , int normalize = 1) 
                tools.AllLegendre(x, maxN, maxN, Plm, dPlm, true, 1, true);
         
                // [Calculate Components of B]

                for (int l = 1; l < maxN + 1; l++) {
                    radialfactor = pow(r, -(l + 2.0));
                    for (int m = 0; m < l + 1; m++) {
                        //-------------find radial b field component---------------
                        gsmdtlm = _COEF.GSMDT[l][m];
                        hsmdtlm = _COEF.HSMDT[l][m];
                        phifactor = gsmdtlm * cos(m * phi) + hsmdtlm * sin(m * phi);
                        Pnm = Plm[l][m];
                        outputBrtp[0] = outputBrtp[0] + ((l + 1) * radialfactor * phifactor * Pnm);

                        //--------------find theta b field component---------------
                        dpnmdtheta = dPlm[l][m]; // dpdtheta(L, M, theta, /inrad)
                        outputBrtp[1] = outputBrtp[1] - radialfactor * phifactor * dpnmdtheta;

                        //-------------------------------------------------------
                        outputBrtp[2] = outputBrtp[2] + m * radialfactor * (gsmdtlm * sin(m * phi) - hsmdtlm * cos(m * phi)) * Pnm;
                    }
                }
                // --------correct bphi and don't divide by 0!------------ 
                sintheta = sin(theta);
                if (sintheta != 0) {
                    outputBrtp[2] /= sintheta;
                } else {
                    outputBrtp[2] = 0.0;
                }
                // [Add CAN model]
                if (CAN) {
                    //  [IAU -> JSM   xyz]
                    vector<DataType> IAUxyz;
                    IAUxyz.resize(3);
                    IAUxyz = tools.Sph2CarP(inputSphCoord_rad);
                    vector<DataType>JSMxyz;
                    JSMxyz.resize(3);
                    JSMxyz = tools.Coord_SIII2JSM(IAUxyz);
                    // [JSM: xyz -> rho z]
                    DataType rho = sqrt(JSMxyz[0] * JSMxyz[0] + JSMxyz[1] * JSMxyz[1]);
                    DataType zzz = JSMxyz[2];
                    vector<DataType> coorIn{rho, zzz};
                    vector<DataType> Bcan; // [rho z, point]
                    Bcan.resize(2);
                    can81(coorIn, Bcan);
                    // [B_JSM: rho z -> xyz]
                    phi = atan2(JSMxyz[1] , JSMxyz[0]);
                    if (JSMxyz[0] < 0)
                        phi += PI;
                    if (phi < 0)
                        phi += 2.0 * PI;
                    vector<DataType> BJSMxyz;
                    BJSMxyz.resize(3);
                    BJSMxyz[0] = Bcan[0] * cos(phi);
                    BJSMxyz[1] = Bcan[0] * sin(phi),
                    BJSMxyz[2] = Bcan[1];
                    // [JSM to System III]
                    vector<DataType> BIAUxyz;
                    BIAUxyz.resize(3);
                    BIAUxyz = tools.Coord_JSM2SIII(BJSMxyz);
                    vector<DataType> BIAUrtp;
                    BIAUrtp.resize(3);
                    BIAUrtp = tools.Car2SphV(IAUxyz, BIAUxyz);
                    outputBrtp[0] = outputBrtp[0] + BIAUrtp[0];
                    outputBrtp[1] = outputBrtp[1] + BIAUrtp[1];
                    outputBrtp[2] = outputBrtp[2] + BIAUrtp[2];
                }
                // [Bx, By, Bz: Bxyz in System III RHS] 
                outputBm = sqrt(outputBrtp[0] * outputBrtp[0] + outputBrtp[1] * outputBrtp[1] + outputBrtp[2] * outputBrtp[2]);
                outputBxyz = tools.Sph2CarV(inputSphCoord_rad, outputBrtp);
                return 0;
            };
           const long double PI = atan(1)*4.0;

        };

        template <typename DataType>
        class Jupiter_JRM09_CAN81BMag : public Jupiter_JRM09_CAN81<DataType, ParamDataSpec<DataType> > {
        public:

            /**
             * @brief Constructor.
             * @details Create the ParamData type of the input ParamData.
             */
            Jupiter_JRM09_CAN81BMag(Process& pProcess, ParamDataSpec<vector<DataType> >& paramInput, bool can) :
            Jupiter_JRM09_CAN81<DataType, ParamDataSpec < DataType > > (pProcess, paramInput, can) {
            }

            virtual ~Jupiter_JRM09_CAN81BMag() {
            }

        protected:

            virtual void pushResult(vector<DataType> /*fieldSph_*/, DataType magnitude_, vector<DataType> /*fieldCart_*/) {
                Jupiter_JRM09_CAN81<DataType, ParamDataSpec<DataType> >::_paramOutput->getDataList().push_back(magnitude_);
            }
        };

        template <typename DataType>
        class Jupiter_JRM09_CAN81Cart : public Jupiter_JRM09_CAN81<DataType, ParamDataSpec<vector<DataType>> >
        {
            public:

            /**
             * @brief Constructor.
             * @details Create the ParamData type of the input ParamData.
             */
            Jupiter_JRM09_CAN81Cart(Process& pProcess, ParamDataSpec<vector<DataType> >& paramInput, bool can) :
                    Jupiter_JRM09_CAN81<DataType, ParamDataSpec <vector< DataType> > > (pProcess, paramInput, can) {
            }

            virtual ~Jupiter_JRM09_CAN81Cart() {
            }

            protected:

            virtual void pushResult(vector<DataType> /*fieldSph_*/, DataType /*magnitude_*/, vector<DataType> fieldCart_) {
                Jupiter_JRM09_CAN81<DataType, ParamDataSpec<vector<DataType> > >::_paramOutput->getDataList().push_back(fieldCart_);
            }
        };

        template <typename DataType>
        class Jupiter_JRM09_CAN81Sphr : public Jupiter_JRM09_CAN81<DataType, ParamDataSpec<vector <DataType> > > {
        public:

            /**
             * @brief Constructor.
             * @details Create the ParamData type of the input ParamData.
             */
            Jupiter_JRM09_CAN81Sphr(Process& pProcess, ParamDataSpec<vector<DataType> >& paramInput, bool can) :
            Jupiter_JRM09_CAN81<DataType, ParamDataSpec <vector <DataType > > > (pProcess, paramInput, can) {
            }

            virtual ~Jupiter_JRM09_CAN81Sphr() {
            }

        protected:

            virtual void pushResult(vector<DataType> fieldSph_, DataType /*magnitude_*/, vector<DataType> /*fieldCart_*/) {
                Jupiter_JRM09_CAN81<DataType, ParamDataSpec<vector<DataType> > >::_paramOutput->getDataList().push_back(fieldSph_);
            }
        };

    } /*en Parameters */
} // end AMDA

#endif /* Jupiter_JRM09_CAN81 */