/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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.
 ******************************************************/

#include "gm_data.h"

gctl::tensor gctl::transform_matrix(const point3ds &op)
{
    tensor R;
    
    R[0][0] = sin((0.5 - op.lat/180.0)*GCTL_Pi)*cos((2.0 + op.lon/180.0)*GCTL_Pi);
    R[0][1] = sin((0.5 - op.lat/180.0)*GCTL_Pi)*sin((2.0 + op.lon/180.0)*GCTL_Pi);
    R[0][2] = cos((0.5 - op.lat/180.0)*GCTL_Pi);

    R[1][0] = cos((0.5 - op.lat/180.0)*GCTL_Pi)*cos((2.0 + op.lon/180.0)*GCTL_Pi);
    R[1][1] = cos((0.5 - op.lat/180.0)*GCTL_Pi)*sin((2.0 + op.lon/180.0)*GCTL_Pi);
    R[1][2] = -1.0*sin((0.5 - op.lat/180.0)*GCTL_Pi);

    R[2][0] = -1.0*sin((2.0 + op.lon/180.0)*GCTL_Pi);
    R[2][1] = cos((2.0 + op.lon/180.0)*GCTL_Pi);
    R[2][2] = 0.0;
    
    return R;
}

void gctl::geomag_local2Cartesian(array<point3dc> &abs_b_, const array<point3dc> &geo_b_, const array<point3ds> &loc_s_)
{
    if (geo_b_.size() != loc_s_.size())
    {
        throw std::runtime_error("[gctl::geomag_local2Cartesian] Incompatible input arraies' size.");
    }
    
    abs_b_.resize(geo_b_.size());

    point3dc mag_v;
    for (size_t i = 0; i < geo_b_.size(); i++)
    {
        // magnetization vector at the local coordinates
		// note here the postive direction of the inclination angle is downward
		// note here the postive direction of the declination angle is clockwise
        mag_v.x = geo_b_[i].z;
		mag_v.y = geo_b_[i].y;
		mag_v.z = -1.0*geo_b_[i].x;

        // magnetic susceptibility is taken as one here
		// rotate the local coordinate system to the regular status
        abs_b_[i] = mag_v.rotate((90.0 - loc_s_[i].lat)*GCTL_Pi/180.0, 0.0, (90.0 + loc_s_[i].lon)*GCTL_Pi/180.0);
    }
    return;
}

void gctl::read_IGRF_table(std::string file, array<IGRF_para> &IGRFs, int head_record, std::string ext_f)
{
    std::ifstream infile;
	gctl::open_infile(infile, file, ext_f);

    std::vector<IGRF_para> vec_para;
    IGRF_para tmp_para;

	std::string err_str, tmp_str;
	std::stringstream tmp_ss;
	std::string date_str, alti_str, dd_str, dm_str, id_str, im_str;

    // 头信息就跳过
    for (size_t i = 0; i < head_record; i++)
    {
        getline(infile, tmp_str);
    }
	
	while (getline(infile, tmp_str))
	{
		//Date Coord-System Altitude Latitude Longitude D_deg D_min I_deg I_min H_nT X_nT Y_nT Z_nT F_nT dD_min dI_min dH_nT dX_nT dY_nT dZ_nT dF_nT
		//2017,6,15 D M4000.00 30.0167 105.017 -2d 16m 46d 52m 34358.9 34332.2 -1354.1 36667.2 50249.6 -3.8 5.6 -13.2 -14.8 -37.9 105.0 67.6
		tmp_ss.clear();
		tmp_ss.str(tmp_str);
		tmp_ss >> date_str >> tmp_para.CSystem >> alti_str 
            >> tmp_para.Latitude >> tmp_para.Longitude
			>> dd_str >> dm_str >> id_str >> im_str
			>> tmp_para.H_nT >> tmp_para.X_nT >> tmp_para.Y_nT 
            >> tmp_para.Z_nT >> tmp_para.F_nT >> tmp_para.dD_min 
            >> tmp_para.dI_min >> tmp_para.dH_nT >> tmp_para.dX_nT 
            >> tmp_para.dY_nT >> tmp_para.dZ_nT >> tmp_para.dF_nT;

		if (1 != sscanf(dd_str.c_str(), "%dd", &tmp_para.D_deg) || 
			1 != sscanf(dm_str.c_str(), "%dm", &tmp_para.D_min) ||
			1 != sscanf(id_str.c_str(), "%dd", &tmp_para.I_deg) ||
			1 != sscanf(im_str.c_str(), "%dm", &tmp_para.I_min) ||
            2 != sscanf(alti_str.c_str(), "%c%lf", &tmp_para.AltiCode, &tmp_para.Altitude))
		{
			throw gctl::invalid_argument("[gctl::read_IGRF_table] Wrong format: " + tmp_str);
		}

        if (3 != sscanf(date_str.c_str(), "%d,%d,%d", &tmp_para.Year, &tmp_para.Month, &tmp_para.Day))
        {
            double deci_year;
            if (1 != sscanf(date_str.c_str(), "%lf", &deci_year))
            {
                throw gctl::invalid_argument("[gctl::read_IGRF_table] Wrong format: " + tmp_str);
            }

            decimal_year2date(deci_year, tmp_para.Year, tmp_para.Month, tmp_para.Day);
        }

		tmp_para.I = fabs(tmp_para.I_deg) + tmp_para.I_min/60.0;
		tmp_para.D = fabs(tmp_para.D_deg) + tmp_para.D_min/60.0;
		if (tmp_para.I_deg < 0) tmp_para.I *= -1;
		if (tmp_para.D_deg < 0) tmp_para.D *= -1;

        vec_para.push_back(tmp_para);
	}
	infile.close();

    IGRFs.input(vec_para);
    destroy_vector(vec_para);
    return;
}

void gctl::magnetic_components2deltaT(const _1d_array &Hax, const _1d_array &Hay, const _1d_array &Za, 
    _1d_array &deltaT, double T0_inclina, double T0_declina)
{
    if (Hax.size() != Hay.size() || Hay.size() != Za.size())
    {
        throw std::runtime_error("[gctl::magnetic_components2deltaT] Unmatched components' size.");
    }
    
    // note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
    double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
    double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
    double za_f  = sin(GCTL_Pi*T0_inclina/180.0);

    deltaT.resize(Hax.size());
    for (size_t i = 0; i < Hax.size(); i++)
    {
        deltaT[i] = hax_f*Hax[i] + hay_f*Hay[i] + za_f*Za[i];
    }
    return;
}

void gctl::magnetic_components2deltaT(const array<point3dc> &Mag_components, _1d_array &deltaT, double T0_inclina, double T0_declina)
{
    // note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
    double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
    double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
    double za_f  = sin(GCTL_Pi*T0_inclina/180.0);

    deltaT.resize(Mag_components.size());
    for (size_t i = 0; i < Mag_components.size(); i++)
    {
        deltaT[i] = hax_f*Mag_components[i].x + hay_f*Mag_components[i].y + za_f*Mag_components[i].z;
    }
    return;
}

void gctl::magnetic_tensors2deltaTs(const array<tensor> &Mag_tensors, array<point3dc> &deltaTs, double T0_inclina, double T0_declina)
{
    // note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
    double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
    double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
    double za_f  = sin(GCTL_Pi*T0_inclina/180.0);

    deltaTs.resize(Mag_tensors.size());
    for (size_t i = 0; i < Mag_tensors.size(); i++)
    {
        deltaTs[i].x = hax_f*Mag_tensors[i][0][0] + hay_f*Mag_tensors[i][0][1] + za_f*Mag_tensors[i][0][2];
        deltaTs[i].y = hax_f*Mag_tensors[i][0][1] + hay_f*Mag_tensors[i][1][1] + za_f*Mag_tensors[i][1][2];
        deltaTs[i].z = hax_f*Mag_tensors[i][0][2] + hay_f*Mag_tensors[i][1][2] + za_f*Mag_tensors[i][2][2];
    }
    return;
}

void gctl::magnetic_components2deltaT_sph(const array<point3dc> &Mag_components, const _1d_array &T0_inclina, const _1d_array &T0_declina, _1d_array &deltaT)
{
    // note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
    deltaT.resize(Mag_components.size());
    for (size_t i = 0; i < Mag_components.size(); i++)
    {
        double hax_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*sin(GCTL_Pi*T0_declina[i]/180.0);
        double hay_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*cos(GCTL_Pi*T0_declina[i]/180.0);
        double za_f  = sin(GCTL_Pi*T0_inclina[i]/180.0);

        // 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)
        deltaT[i] = hax_f*Mag_components[i].z - hay_f*Mag_components[i].y + za_f*Mag_components[i].x;
    }
    return;
}

double gctl::magnetic_components2deltaT_sph(const point3dc &Mag_components, double T0_inclina, double T0_declina)
{
    double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
    double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
    double za_f  = sin(GCTL_Pi*T0_inclina/180.0);

    // 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)
    return hax_f*Mag_components.z - hay_f*Mag_components.y + za_f*Mag_components.x;
}

void gctl::magnetic_tensors2deltaTs_sph(const array<tensor> &Mag_tensors, const _1d_array &T0_inclina, const _1d_array &T0_declina, array<point3dc> &deltaTs)
{
    // note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
    deltaTs.resize(Mag_tensors.size());
    for (size_t i = 0; i < Mag_tensors.size(); i++)
    {
        double hax_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*sin(GCTL_Pi*T0_declina[i]/180.0);
        double hay_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*cos(GCTL_Pi*T0_declina[i]/180.0);
        double za_f  = sin(GCTL_Pi*T0_inclina[i]/180.0);

        // note here Mag_tensors[i][:][0] is the reversed radial component
        // Mag_components[i][:][1] is the latitudinal component (south pointing). to use it in a local cratesian coordinate, we need to reverse it
        // Mag_components[i][:][2] is the longtidinal component (east pointing)
        deltaTs[i].x = za_f*Mag_tensors[i][0][0] - hay_f*Mag_tensors[i][0][1] + hax_f*Mag_tensors[i][0][2];
        deltaTs[i].y = za_f*Mag_tensors[i][0][1] - hay_f*Mag_tensors[i][1][1] + hax_f*Mag_tensors[i][1][2];
        deltaTs[i].z = za_f*Mag_tensors[i][0][2] - hay_f*Mag_tensors[i][1][2] + hax_f*Mag_tensors[i][2][2];
    }
    return;
}

void gctl::power_spectrum_shc(array<double> &P, const array<double> &C, const array<double> &S, double R, double r, int n, int N)
{
    if (n < 0 || n >= N || C.size() < (N + 1)*(N + 2)/2 || S.size() < (N + 1)*(N + 2)/2)
    {
        throw std::runtime_error("[gctl::power_spectrum_shc] Invalid input parameters.");
    }

    P.resize(N - n + 1, 0);

    int id;
    for (int i = n; i < N + 1; i++)
    {
        id = i*(i + 1)/2;
        for (int j = 0; j < i + 1; j++)
        {
            P[i - n] += (C[id + j]*C[id + j] + S[id + j]*S[id + j]);
        }

        P[i - n] *= (n + 1)*pow(R/r, 2*n + 4);
    }
    return;
}