/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_GM_DATA_H
#define _GCTL_GM_DATA_H

#include "gctl_potential_config.h"
#ifdef GCTL_POTENTIAL_OPENMP
#include "omp.h"
#include <thread>
#endif

#include "gctl/core.h"
#include "gctl/utility.h"
#include "gctl/geometry.h"
#include "gctl/maths.h"
#include "gctl/algorithms.h"

namespace gctl
{
    /**
     * @brief      components of potential field
     */
    enum gravitational_field_type_e
    {
        GravPot, ///< field potential
        Vx, ///< x-gradient of V (longitudinal-gradient in the spherical coordinates)
        Vy, ///< y-gradient of V (latitudinal-gradient in the spherical coordinates)
        Vz, ///< z-gradient of V (radial-gradient in the spherical coordinates)
        Txx, ///< x-gradient of Gx
        Txy, ///< y-gradient of Gx
        Txz, ///< z-gradient of Gx
        Tyx, ///< x-gradient of Gy
        Tyy, ///< y-gradient of Gy
        Tyz, ///< z-gradient of Gy
        Tzx, ///< x-gradient of Gz
        Tzy, ///< y-gradient of Gz
        Tzz, ///< z-gradient of Gz
    };

    enum magnetic_field_type_e
    {
        MagPot,
        Hax,
        Hay,
        Za,
        Bx,
        By,
        Bz,
        Bxx,
        Bxy,
        Bxz,
        Byx,
        Byy,
        Byz,
        Bzx,
        Bzy,
        Bzz,
        DeltaT, ///< total magnetic strength
        DeltaTx, ///< x-gradient of total magnetic strength
        DeltaTy, ///< x-gradient of total magnetic strength
        DeltaTz, ///< x-gradient of total magnetic strength
    };

    /**
     * @brief A structure holds all outputs of the standard IGRF or EMM program.
     * 
     */
    struct IGRF_para
    {
        int Year, Month, Day;
        char CSystem, AltiCode;
        double Altitude, Latitude, Longitude;
        double D, I, H_nT, X_nT, Y_nT, Z_nT, F_nT, dD_min, dI_min, dH_nT, dX_nT, dY_nT, dZ_nT, dF_nT; 
        int D_deg, D_min, I_deg, I_min;
    };

    /**
     * @brief Calculate the transform matrix of the vector componments wrt. the spherical coordinates.
     * 
     * @param op The observation point on a sphere
     * @return The transform matrix
     */
    tensor transform_matrix(const point3ds &op);

    /**
     * @brief Transform localized magnetization vectors on a sphere to the absolute Cartesian coordinates.
     * 
     * @param abs_b_ Output absolute vectors
     * @param geo_b_ Local magnetization vectors on a sphere
     * @param loc_s_ Spherical coordinates of the local vectors
     */
    void geomag_local2Cartesian(array<point3dc> &abs_b_, const array<point3dc> &geo_b_, const array<point3ds> &loc_s_);

    /**
     * @brief Read output table of the standard IGRF or EMM program
     * 
     * @param file Input filename
     * @param IGRFs Return IGRF parameters
     * @param head_record Lines of head records
     * @param ext_f File name extension
     */
    void read_IGRF_table(std::string file, array<IGRF_para> &IGRFs, int head_record = 1, std::string ext_f = ".txt");

    /**
     * @brief Transform magnetic components data to delta_T anomalies data.
     * 
     * @param Hax Hax component of the magnetic data
     * @param Hay Hay component of the magnetic data
     * @param Za Za component of the magnetic data
     * @param deltaT Output delta_T anomalies data
     * @param T0_inclina Inclination degree of the earth normal magnetic field.
     * @param T0_declina Declination degree of the earth normal magnetic field.
     */
    void magnetic_components2deltaT(const _1d_array &Hax, const _1d_array &Hay, const _1d_array &Za, 
        _1d_array &deltaT, double T0_inclina, double T0_declina);

    /**
     * @brief Transform magnetic components data to delta_T anomalies data.
     * 
     * @param Mag_components Magnetic components of the magnetic data
     * @param deltaT Output delta_T anomalies data
     * @param T0_inclina Inclination degree of the earth normal magnetic field.
     * @param T0_declina Declination degree of the earth normal magnetic field.
     */
    void magnetic_components2deltaT(const array<point3dc> &Mag_components, _1d_array &deltaT, double T0_inclina, double T0_declina);

    /**
     * @brief Transform magnetic tensor data to delta_T gradient anomalies data.
     * 
     * @param Mag_tensors Magnetic tensors of the magnetic data
     * @param deltaTs Output delta_T gradient anomalies data
     * @param T0_inclina Inclination degree of the earth normal magnetic field.
     * @param T0_declina Declination degree of the earth normal magnetic field.
     */
    void magnetic_tensors2deltaTs(const array<tensor> &Mag_tensors, array<point3dc> &deltaTs, double T0_inclina, double T0_declina);

    /**
     * @brief Transform magnetic components data to delta_T anomalies data.
     * 
     * @note note here Mag_components[i].x is the reversed radial component. 
     * Mag_components[i].y is the latitudinal component (south pointing). 
     * to use it in a local cratesian coordinate, we need to reverse it. 
     * Mag_components[i].z is the longtidinal component (east pointing)
     * 
     * @param Mag_components Magnetic components of the magnetic data
     * @param deltaT Output delta_T anomalies data
     * @param T0_inclina Inclination degrees of the earth normal magnetic field.
     * @param T0_declina Declination degrees of the earth normal magnetic field.
     */
    void magnetic_components2deltaT_sph(const array<point3dc> &Mag_components, const _1d_array &T0_inclina, const _1d_array &T0_declina, _1d_array &deltaT);

    /**
     * @brief Transform magnetic components data to delta_T anomalies data.
     * 
     * @note note here Mag_components.x is the reversed radial component. 
     * Mag_components.y is the latitudinal component (south pointing). 
     * to use it in a local cratesian coordinate, we need to reverse it. 
     * Mag_components.z is the longtidinal component (east pointing)
     * 
     * @param Mag_components Magnetic components of the magnetic data
     * @param T0_inclina Inclination degrees of the earth normal magnetic field.
     * @param T0_declina Declination degrees of the earth normal magnetic field.
     */
    double magnetic_components2deltaT_sph(const point3dc &Mag_components, double T0_inclina, double T0_declina);

    /**
     * @brief Transform magnetic tensor data to delta_T gradient anomalies data.
     * 
     * @param Mag_tensors Magnetic tensors of the magnetic data
     * @param deltaTs Output delta_T gradient anomalies data
     * @param T0_inclina Inclination degrees of the earth normal magnetic field.
     * @param T0_declina Declination degrees of the earth normal magnetic field.
     */
    void magnetic_tensors2deltaTs_sph(const array<tensor> &Mag_tensors, const _1d_array &T0_inclina, const _1d_array &T0_declina, array<point3dc> &deltaTs);

    /**
     * @brief Calculate the power spectrum of given spherical harmonic coefficients.
     * 
     * @param P Returned power spectrum
     * @param C Cosine coefficients strats from 0 order 0 degree up to least N order N degree
     * @param S Sine coefficients strats from 0 order 0 degree up to least N order N degree
     * @param R Reference radius (usually is 6371.2 km)
     * @param r Observation radius
     * @param n Strat order
     * @param N End order
     */
    void power_spectrum_shc(array<double> &P, const array<double> &C, const array<double> &S, 
        double R, double r, int n, int N);
}

#endif // _GCTL_GM_DATA_H