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

void gctl::mag_axial_dipole_parameters(double &M, double &I, double r, double lat_deg, double sus, double ref_r, double abs_g0)
{
	M = sus*fabs(abs_g0)*power3(ref_r/r)*sqrt(1.0 + 3.0*power2(sind(lat_deg)));
	I = atan(2.0*tan(GCTL_Pi*lat_deg/180.0));
	return;
}

gctl::point3dc gctl::magkernel_single(const mag_dipole &a_dipole, const point3dc &a_op)
{
    double r = a_op.module();
    // T -> nT 1e+9 * 1e-7 -> 1e+2
    return 1e+2*a_dipole.M*(3.0*dot(a_dipole.n, a_op)/power5(r)*a_op - 1.0/power3(r)*a_dipole.n);
}

gctl::point3dc gctl::magkernel_single(const mag_dipole &a_dipole, const point3ds &a_op, tensor *R_ptr)
{
    tensor R;
    if (R_ptr != nullptr) R = *R_ptr;
	else R = transform_matrix(a_op);

    double r = a_op.rad;
    point3dc pc = a_op.s2c();

    // T -> nT 1e+9 * 1e-7 -> 1e+2
    point3dc b = 1e+2*a_dipole.M*(3.0*dot(a_dipole.n, pc)/power5(r)*pc - 1.0/power3(r)*a_dipole.n);
    b = R*b;
	b.z = 0.0; // the latitudinal componement is zero
    return b;
}

void gctl::magobser(array<point3dc> &out_obs, const array<mag_dipole> &dipoles, 
    const array<point3dc> &obsp, verbose_type_e verbose)
{
    int i, j;
	int o_size = obsp.size();
	int e_size = dipoles.size();

	out_obs.resize(o_size, point3dc(0.0, 0.0, 0.0));

	gctl::progress_bar bar(e_size, "magobser_componments");
	for (j = 0; j < e_size; j++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(j);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);

#pragma omp parallel for private (i) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			out_obs[i] = out_obs[i] + magkernel_single(dipoles[j], obsp[i]);;
		}
	}
    return;
}

void gctl::magobser(array<point3dc> &out_obs, const array<mag_dipole> &dipoles, 
    const array<point3ds> &obsp, verbose_type_e verbose)
{
    int i, j;
	int o_size = obsp.size();
	int e_size = dipoles.size();

	out_obs.resize(o_size, point3dc(0.0, 0.0, 0.0));
	array<tensor> Rs(o_size);
#pragma omp parallel for private (i) schedule(guided)
	for (i = 0; i < o_size; i++)
	{
		Rs[i] = transform_matrix(obsp[i]);
	}

	gctl::progress_bar bar(e_size, "magobser_componments");
	for (j = 0; j < e_size; j++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(j);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);

#pragma omp parallel for private (i) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			out_obs[i] = out_obs[i] + magkernel_single(dipoles[j], obsp[i], Rs.get(i));
		}
	}
    return;
}