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

void gctl::magobser(array<point3dc> &out_obs, const shc_data &sd, const array<point3ds>& obsp, verbose_type_e verbose)
{
    int id;
    double rad, fx, fy, fz;
    _1d_array P, dP;

    out_obs.resize(obsp.size(), point3dc(0.0, 0.0, 0.0));

    progress_bar bar(obsp.size(), "magobser_shc");
    for (int p = 0; p < obsp.size(); p++)
    {
        if (verbose == gctl::FullMsg) bar.progressed(p);
        else if (verbose == gctl::ShortMsg) bar.progressed_simple(p);

        rad = obsp[p].rad;
        schmidt_legendre(90.0 - obsp[p].lat, 450, P, dP);

        // 注意跳过第0阶0次项
        for (int i = sd.ns; i <= sd.ne; i++)
        {
            id = i*(i + 1)/2;
            fx = pow(sd.rad/rad, i+2);
            fz = -1.0*(i+1)*pow(sd.rad/rad, i+2);

            for (int j = 0; j <= i; j++)
            {
                fy = j/sind(90.0 - obsp[p].lat)*pow(sd.rad/rad, i+2);

                out_obs[p].x += -1.0*fx*(cosd(j*obsp[p].lon)*dP[id + j]*sd.coeff_c(i, j) + sind(j*obsp[p].lon)*dP[id + j]*sd.coeff_s(i, j));
                out_obs[p].y += fy*(sind(j*obsp[p].lon)*P[id + j]*sd.coeff_c(i, j) - cosd(j*obsp[p].lon)*P[id + j]*sd.coeff_s(i, j));
                out_obs[p].z += fz*(cosd(j*obsp[p].lon)*P[id + j]*sd.coeff_c(i, j) + sind(j*obsp[p].lon)*P[id + j]*sd.coeff_s(i, j));
            }
        }
    }
    return;
}