/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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 = abs(tmp_para.I_deg) + tmp_para.I_min/60.0;
		tmp_para.D = abs(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::read_Swarm_shc(std::string file, array<shc_data> &SHCs, int &spline_odr, int &sample_num)
{
    dsv_io tio;
    tio.head_number(2);
    tio.load_text(file);

    // parse the head records
    int n, N, tnum;
    std::vector<double> time_snaps;
    std::vector<std::string> rcd = tio.head_records();
    parse_string_to_value(rcd[0], ' ', true, n, N, tnum, spline_odr, sample_num);
    parse_string_to_vector(rcd[1], ' ', time_snaps);

    array<int> nid, mid;
    array<double> vals;
    int rnum = tio.row_number();
    tio.get_column(nid, 1);
    tio.get_column(mid, 2);

    // store coefficients
    int id, shc_num;
    SHCs.resize(tnum);
    for (size_t i = 0; i < tnum; i++)
    {
        SHCs[i].ns = n;
        SHCs[i].ms = 0;
        SHCs[i].ne = N;
        SHCs[i].me = N;
        SHCs[i].time = time_snaps[i];

        shc_num = (N + 1)*(N + 2)/2;
        SHCs[i].Snm.resize(shc_num, 0.0);
        SHCs[i].Cnm.resize(shc_num, 0.0);

        tio.get_column(vals, i + 3);
        for (size_t j = 0; j < rnum; j++)
        {
            if (mid[j] > 0)
            {
                id = nid[j]*(nid[j] + 1)/2 + mid[j];
                SHCs[i].Snm[id] = vals[j];
                SHCs[i].Cnm[id] = vals[j + 1];
            }
            else if (mid[j] == 0)
            {
                id = nid[j]*(nid[j] + 1)/2;
                SHCs[i].Snm[id] = vals[j];
                SHCs[i].Cnm[id] = 0.0;
            }
        }
    }
    return;
}

void gctl::save_Swarm_shc(std::string file, const array<shc_data> &SHCs, int spline_odr, int sample_num, std::string info)
{
    time_t now = time(0);
	char* dt = ctime(&now);

    std::ofstream ofile;
    open_outfile(ofile, file + ".shc", ".txt");
    ofile << "# " << info << std::endl;
    ofile << "# Generated on " << dt;

    int s = SHCs[0].ns, S = SHCs[0].ne;
    int c = SHCs.size();
    ofile << s << " " << S << " " << c << " " << spline_odr << " " << sample_num << std::endl;

    ofile << SHCs[0].time;
    for (size_t i = 1; i < SHCs.size(); i++)
    {
        ofile << " " << SHCs[i].time;
    }

    for (size_t n = s; n <= S; n++)
    {
        for (size_t m = 0; m <= n; m++)
        {
            if (m == 0)
            {
                ofile << "\n" << n << " " << m;
                for (size_t i = 0; i < c; i++)
                {
                    ofile << " " << SHCs[i].coeff_s(n, m);
                }
            }
            else
            {
                ofile << "\n" << n << " " << m;
                for (size_t i = 0; i < c; i++)
                {
                    ofile << " " << SHCs[i].coeff_s(n, m);
                }

                ofile << "\n" << n << " -" << m;
                for (size_t i = 0; i < c; i++)
                {
                    ofile << " " << SHCs[i].coeff_c(n, m);
                }
            }
        }
    }
    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::schmidt_legendre(double co_lati, int n_max, array<double> &P, array<double> &dP)
{
    if (n_max <= 0) throw std::runtime_error("[gctl::schmidt_legendre] Invalid paramters.");

    double cos_lat = cosd(co_lati);
    if (abs(cos_lat) == 1.0) throw std::runtime_error("[gctl::schmidt_legendre] Derivative cannot be calculated at poles.");

    double pm2, pm1, pmm, plm, rescalem, z, scalef;
    int k, kstart, m, n, NumTerms;

    NumTerms = ((n_max + 1) * (n_max + 2) / 2);
    P.resize(NumTerms, 0.0);
    dP.resize(NumTerms, 0.0);
    _1d_array f1(NumTerms + 1), f2(NumTerms + 1), PreSqr(NumTerms + 1);

    scalef = 1.0e-280;

    for(n = 0; n <= 2 * n_max + 1; ++n)
    {
        PreSqr[n] = sqrt((double) (n));
    }

    k = 2;
    for(n = 2; n <= n_max; n++)
    {
        k = k + 1;
        f1[k] = (double) (2 * n - 1) / (double) (n);
        f2[k] = (double) (n - 1) / (double) (n);
        for(m = 1; m <= n - 2; m++)
        {
            k = k + 1;
            f1[k] = (double) (2 * n - 1) / PreSqr[n + m] / PreSqr[n - m];
            f2[k] = PreSqr[n - m - 1] * PreSqr[n + m - 1] / PreSqr[n + m] / PreSqr[n - m];
        }
        k = k + 2;
    }

    /*z = sin (geocentric latitude) */
    z = sqrt((1.0 - cos_lat)*(1.0 + cos_lat));
    pm2 = 1.0;
    P[0] = 1.0;
    dP[0] = 0.0;

    pm1 = cos_lat;
    P[1] = pm1;
    dP[1] = z;
    k = 1;

    for(n = 2; n <= n_max; n++)
    {
        k = k + n;
        plm = f1[k] * cos_lat * pm1 - f2[k] * pm2;
        P[k] = plm;
        dP[k] = (double) (n) * (pm1 - cos_lat * plm) / z;
        pm2 = pm1;
        pm1 = plm;
    }

    pmm = PreSqr[2] * scalef;
    rescalem = 1.0 / scalef;
    kstart = 0;

    for(m = 1; m <= n_max - 1; ++m)
    {
        rescalem = rescalem*z;

        /* Calculate Pcup(m,m)*/
        kstart = kstart + m + 1;
        pmm = pmm * PreSqr[2 * m + 1] / PreSqr[2 * m];
        P[kstart] = pmm * rescalem / PreSqr[2 * m + 1];
        dP[kstart] = -((double) (m) * cos_lat * P[kstart] / z);
        pm2 = pmm / PreSqr[2 * m + 1];
        /* Calculate Pcup(m+1,m)*/
        k = kstart + m + 1;
        pm1 = cos_lat * PreSqr[2 * m + 1] * pm2;
        P[k] = pm1*rescalem;
        dP[k] = ((pm2 * rescalem) * PreSqr[2 * m + 1] - cos_lat * (double) (m + 1) * P[k]) / z;
        /* Calculate Pcup(n,m)*/
        for(n = m + 2; n <= n_max; ++n)
        {
            k = k + n;
            plm = cos_lat * f1[k] * pm1 - f2[k] * pm2;
            P[k] = plm*rescalem;
            dP[k] = (PreSqr[n + m] * PreSqr[n - m] * (pm1 * rescalem) - (double) (n) * cos_lat * P[k]) / z;
            pm2 = pm1;
            pm1 = plm;
        }
    }

    /* Calculate Pcup(nMax,nMax)*/
    rescalem = rescalem*z;
    kstart = kstart + m + 1;
    pmm = pmm / PreSqr[2 * n_max];
    P[kstart] = pmm * rescalem;
    dP[kstart] = -(double) (n_max) * cos_lat * P[kstart] / z;
    return;
}

void gctl::power_spectrum(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] 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] *= (i + 1)*pow(R/r, 2*i + 4);
    }
    return;
}