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


/**
 * @brief      callback interface of the gravitational kernel of 2D triangular elements
 *
 * @param      out_kernel  The returned gravitational kernel
 * @param[in]  ele         2D triangular elements
 * @param[in]  obsp        2D observation points
 * @param[in]  verbose     The verbose level
 */
typedef void (*gkernel_triangle2d_ptr)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);

void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vxz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose);

/**
 * @brief      Integrated callback function of the gravitational kernel of 2D triangle elements
 *
 * @param      out_kernel  The returned gravitational kernel
 * @param[in]  ele         2D triangular elements
 * @param[in]  obsp        2D observation points
 * @param[in]  comp_id     The gravitational component identifier
 * @param[in]  verbose     The verbose level
 */
void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele, 
	const array<point2dc> &obsp, gravitational_field_type_e comp_id, verbose_type_e verbose)
{
	gkernel_triangle2d_ptr triangle_kernel;
	switch (comp_id)
	{
		case Vz:
			triangle_kernel = gkernel_triangle2d_vz;
			break;
		case Vx:
			triangle_kernel = gkernel_triangle2d_vx;
			break;
		case Txx:
			triangle_kernel = gkernel_triangle2d_vxx;
			break;
		case Txz:
			triangle_kernel = gkernel_triangle2d_vxz;
			break;
		case Tzx:
			triangle_kernel = gkernel_triangle2d_vzx;
			break;
		case Tzz:
			triangle_kernel = gkernel_triangle2d_vzz;
			break;
		default:
			throw std::invalid_argument("Invalid gravitational field type!");
			break;
	}
	return triangle_kernel(out_kernel, ele, obsp, verbose);
}

typedef void (*gkernel_triangle2d_ptr2)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);

void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vxx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vxz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);
void gkernel_triangle2d_vzz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose);

void gctl::gkernel(matrix<double> &out_kernel, const array<triangle2d> &ele, 
	const array<point2dp> &obsp, gravitational_field_type_e comp_id, verbose_type_e verbose)
{
	gkernel_triangle2d_ptr2 triangle_kernel2;
	switch (comp_id)
	{
		case Vz:
			triangle_kernel2 = gkernel_triangle2d_vz2;
			break;
		case Vx:
			triangle_kernel2 = gkernel_triangle2d_vx2;
			break;
		case Txx:
			triangle_kernel2 = gkernel_triangle2d_vxx2;
			break;
		case Txz:
			triangle_kernel2 = gkernel_triangle2d_vxz2;
			break;
		case Tzx:
			triangle_kernel2 = gkernel_triangle2d_vzx2;
			break;
		case Tzz:
			triangle_kernel2 = gkernel_triangle2d_vzz2;
			break;
		default:
			throw std::invalid_argument("Invalid gravitational field type!");
			break;
	}
	return triangle_kernel2(out_kernel, ele, obsp, verbose);
}

/**
 * @brief      callback interface of the gravitational observations of 2D triangular elements
 *
 * @param      out_obs     The returned gravitational data
 * @param[in]  ele         2D triangular elements
 * @param[in]  obsp        2D observation points
 * @param[in]  rho         Densities of the triangular elements
 * @param[in]  verbose     The verbose level
 */
typedef void (*gobser_triangle2d_ptr)(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);

void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vxz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);

/**
 * @brief      Integrated callback function of the gravitational observations of 2D triangle elements
 *
 * @param      out_kernel  The returned gravitational kernel
 * @param[in]  ele         2D triangular elements
 * @param[in]  obsp        2D observation points
 * @param[in]  rho         Densities of the triangular elements
 * @param[in]  comp_id     The gravitational component identifier
 * @param[in]  verbose     The verbose level
 */
void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dc> &obsp, 
	const array<double> &rho, gravitational_field_type_e comp_id, verbose_type_e verbose)
{
	gobser_triangle2d_ptr triangle_obser;
	switch (comp_id)
	{
		case Vz:
			triangle_obser = gobser_triangle2d_vz;
			break;
		case Vx:
			triangle_obser = gobser_triangle2d_vx;
			break;
		case Txx:
			triangle_obser = gobser_triangle2d_vxx;
			break;
		case Txz:
			triangle_obser = gobser_triangle2d_vxz;
			break;
		case Tzx:
			triangle_obser = gobser_triangle2d_vzx;
			break;
		case Tzz:
			triangle_obser = gobser_triangle2d_vzz;
			break;
		default:
			throw std::invalid_argument("Invalid gravitational field type!");
			break;
	}
	return triangle_obser(out_obs, ele, obsp, rho, verbose);
}

typedef void (*gobser_triangle2d_ptr2)(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);

void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vxx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vxz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);
void gobser_triangle2d_vzz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose);

void gctl::gobser(array<double> &out_obs, const array<triangle2d> &ele, const array<point2dp> &obsp, 
	const array<double> &rho, gravitational_field_type_e comp_id, verbose_type_e verbose)
{
	gobser_triangle2d_ptr2 triangle_obser2;
	switch (comp_id)
	{
		case Vz:
			triangle_obser2 = gobser_triangle2d_vz2;
			break;
		case Vx:
			triangle_obser2 = gobser_triangle2d_vx2;
			break;
		case Txx:
			triangle_obser2 = gobser_triangle2d_vxx2;
			break;
		case Txz:
			triangle_obser2 = gobser_triangle2d_vxz2;
			break;
		case Tzx:
			triangle_obser2 = gobser_triangle2d_vzx2;
			break;
		case Tzz:
			triangle_obser2 = gobser_triangle2d_vzz2;
			break;
		default:
			throw std::invalid_argument("Invalid gravitational field type!");
			break;
	}
	return triangle_obser2(out_obs, ele, obsp, rho, verbose);
}

double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);
double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr);

/** 我们在这里使用的坐标系（与建模环境兼容）
 * y
 * |
 * |
 * |
 * --------------> x
 * 
 * 正演公式使用的坐标系
 * --------------> x
 * |
 * |
 * |
 * z
 * 
 * z和y含义相同（第二正交轴）但方向相反，上下两个坐标系互成镜像。我们使用的三角形排列为逆时针，相当于正演坐标系中的顺时针，所以得到值为负。
 * 需要再乘以负号反转。
 * 
 * A=((x[j]-i)*z[j+1]-(x[j+1]-i)*z[j])/(pow(z[j+1]-z[j],2)+pow(x[j+1]-x[j],2));	
 * B=(x[j+1]-x[j])*(arctg(z[j]/(x[j]-i),x[j+1]-x[j])-arctg(z[j+1]/(x[j+1]-i),x[j+1]-x[j]));
 * C=0.5*(z[j+1]-z[j])*log((pow((x[j+1]-i),2)+z[j+1]*z[j+1])/(pow((x[j]-i),2)+z[j]*z[j]));
 * g+=A*(B+C);
 * g*=2000*G*den;
 */
void gkernel_triangle2d_vz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vxx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vxx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vxz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vxz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vzx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vzx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vzz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vzz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

#pragma omp parallel for private (j) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			out_kernel[i][j] = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i));
		}
	}
	return;
}

void gkernel_triangle2d_vz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc, obsg;
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, obsg) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
			obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
			out_kernel[i][j] = obsg.y*sin(obsp[i].arc) - obsg.x*cos(obsp[i].arc);
		}
	}
}

void gkernel_triangle2d_vx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc, obsg;
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, obsg) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc);
			obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc);
			out_kernel[i][j] = -1.0*obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc);
		}
	}
}

void gkernel_triangle2d_vxx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vxx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, t) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_kernel[i][j] = sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc)) 
				- cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc));
		}
	}
}

void gkernel_triangle2d_vxz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vxz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, t) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_kernel[i][j] = -1.0*cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc)) 
				+ sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc));
		}
	}
}

void gkernel_triangle2d_vzx2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vzx");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, t) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_kernel[i][j] = sin(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc)) 
				- cos(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc));
		}
	}
}

void gkernel_triangle2d_vzz2(gctl::matrix<double> &out_kernel, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_kernel.resize(o_size, e_size);

	gctl::progress_bar bar(o_size, "gkernel_vzz");
	for (i = 0; i < o_size; i++)
	{
		if (verbose == gctl::FullMsg) bar.progressed(i);
		else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);

		obsc = obsp[i].p2c();

#pragma omp parallel for private (j, t) schedule(guided)
		for (j = 0; j < e_size; j++)
		{
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_kernel[i][j] = -1.0*cos(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc)) 
				+ sin(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc));
		}
	}
}

void gobser_triangle2d_vz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vz");
	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] += gkernel_triangle2d_vz_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vx");
	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] += gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vxx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vxx");
	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] += gkernel_triangle2d_vxx_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vxz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vxz");
	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] += gkernel_triangle2d_vxz_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vzx");
	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] += gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dc> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vzz");
	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] += gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)) * rho[j];
		}
	}
	return;
}

void gobser_triangle2d_vz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc, obsg;
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vz");
	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, obsg, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc); // 沿x轴正方向的重力值
			obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc); // 沿z轴负方向的重力值
			out_obs[i] += (obsg.y*sin(obsp[i].arc) - obsg.x*cos(obsp[i].arc))*rho[j]; // 沿半径负方向的重力值
		}
	}
	return;
}

void gobser_triangle2d_vx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc, obsg;
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vx");
	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, obsg, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			obsg.x = gkernel_triangle2d_vx_sig(ele.get(j), &obsc); // 沿x轴正方向的重力值
			obsg.y = gkernel_triangle2d_vz_sig(ele.get(j), &obsc); // 沿z轴负方向的重力值
			out_obs[i] += (-1.0*obsg.y*cos(obsp[i].arc) + obsg.x*sin(obsp[i].arc))*rho[j]; // 沿弧度负方向的重力值
		}
	}
	return;
}

void gobser_triangle2d_vxx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vxx");
	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, t, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_obs[i] += (sin(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc)) 
				- cos(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc)))*rho[j];
		}
	}
}

void gobser_triangle2d_vxz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vxz");
	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, t, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_obs[i] += (-1.0*cos(obsp[i].arc)*(t[0][0]*sin(obsp[i].arc) - t[1][0]*cos(obsp[i].arc)) 
				+ sin(obsp[i].arc)*(t[0][1]*sin(obsp[i].arc) - t[1][1]*cos(obsp[i].arc)))*rho[j];
		}
	}
}

void gobser_triangle2d_vzx2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vzx");
	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, t, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_obs[i] += (sin(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc)) 
				- cos(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc)))*rho[j];
		}
	}
}

void gobser_triangle2d_vzz2(gctl::array<double> &out_obs, const gctl::array<gctl::triangle2d> &ele, 
	const gctl::array<gctl::point2dp> &obsp, const gctl::array<double> &rho, gctl::verbose_type_e verbose)
{
	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

	gctl::point2dc obsc;
	double t[2][2];
	out_obs.resize(o_size, 0.0);

	gctl::progress_bar bar(e_size, "gobser_vzz");
	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, t, obsc) schedule(guided)
		for (i = 0; i < o_size; i++)
		{
			obsc = obsp[i].p2c();
			t[1][0] = gkernel_triangle2d_vzx_sig(ele.get(j), &obsc);
			t[1][1] = gkernel_triangle2d_vzz_sig(ele.get(j), &obsc);
			t[0][0] = -1.0*t[1][1];
			t[0][1] = t[1][0];
			out_obs[i] += (-1.0*cos(obsp[i].arc)*(t[1][0]*sin(obsp[i].arc) - t[0][0]*cos(obsp[i].arc)) 
				+ sin(obsp[i].arc)*(t[1][1]*sin(obsp[i].arc) - t[0][1]*cos(obsp[i].arc)))*rho[j];
		}
	}
}

double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	double sum;
	double A, B, C;
	double Aa, Ab, Ba, Bb, Ca, Cb;

	sum = 0.0;
	for (int n = 0; n < 3; n++)
	{
		Aa = (ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y) -
			 (ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[n]->y - op_ptr->y);
		Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
			 pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);
		A = Aa/Ab;

		Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
		Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
		if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
		if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
		B = (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba - Bb);

		Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
		Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
		C = 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(log(Ca) - log(Cb));

		sum += A*(B + C);
	}
	return -2.0e+8*GCTL_G0*sum;
}

double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	double sum;
	double A, B, C;
	double Aa, Ab, Ba, Bb, Ca, Cb;

	sum = 0.0;
	for (int n = 0; n < 3; n++)
	{
		Aa = (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) -
			 (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[n]->x - op_ptr->x);
		Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
			 pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);
		A = Aa/Ab;

		Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
		Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
		if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
		if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
		B = (ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ba - Bb);

		Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
		Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
		C = 0.5*(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(log(Cb) - log(Ca));

		sum += A*(B + C);
	}
	return 2.0e+8*GCTL_G0*sum;
}

double gkernel_triangle2d_vxx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	return -1.0*gkernel_triangle2d_vzz_sig(ele_ptr, op_ptr);
}

double gkernel_triangle2d_vxz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	return gkernel_triangle2d_vzx_sig(ele_ptr, op_ptr);
}

double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	double sum;
	double A, Ax, B, Bx, C, Cx;
	double Aa, Aa_x, Ab, Ba, Ba_x, Bb, Bb_x, Ca, Ca_x, Cb, Cb_x;
	double tmp_d;

	sum = 0.0;
	for (int n = 0; n < 3; n++)
	{
		Aa = (ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y) -
			 (ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[n]->y - op_ptr->y);
		Aa_x = ele_ptr->vert[n]->y - ele_ptr->vert[(n+1)%3]->y;
		Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
			 pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);

		A = Aa/Ab;
		Ax= Aa_x/Ab;

		Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
		if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
		Ba_x = (ele_ptr->vert[n]->y - op_ptr->y)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x) 
			+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));

		Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
		if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
		Bb_x = (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) 
			+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));

		B = (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba - Bb);
		Bx= (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba_x - Bb_x);

		Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
		Ca_x = -2.0*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
		Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
		Cb_x = -2.0*(ele_ptr->vert[n]->x - op_ptr->x);

		C = 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(log(Ca) - log(Cb));
		Cx= 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ca_x/Ca - Cb_x/Cb);

		sum += (Ax*(B+C) + A*(Bx+Cx));
	}
	return -2.0e+8*GCTL_G0*sum;
}

double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr)
{
	double sum;
	double A, Ay, B, By, C, Cy;
	double Aa, Aa_y, Ab, Ba, Ba_y, Bb, Bb_y, Ca, Ca_y, Cb, Cb_y;
	double tmp_d;

	sum = 0.0;
	for (int n = 0; n < 3; n++)
	{
		Aa = (ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y) -
			 (ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[n]->y - op_ptr->y);
		Aa_y = ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x;
		Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) +
			 pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2);

		A = Aa/Ab;
		Ay= Aa_y/Ab;

		Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x);
		if (ele_ptr->vert[n]->y - op_ptr->y > 0.0 && ele_ptr->vert[n]->x - op_ptr->x < 0.0) Ba -= 2*GCTL_Pi;
		Ba_y = -1.0*(ele_ptr->vert[n]->x - op_ptr->x)/((ele_ptr->vert[n]->x - op_ptr->x)*(ele_ptr->vert[n]->x - op_ptr->x) 
			+ (ele_ptr->vert[n]->y - op_ptr->y)*(ele_ptr->vert[n]->y - op_ptr->y));

		Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x);
		if (ele_ptr->vert[(n+1)%3]->y - op_ptr->y > 0.0 && ele_ptr->vert[(n+1)%3]->x - op_ptr->x < 0.0) Bb -= 2*GCTL_Pi;
		Bb_y = -1.0*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x)/((ele_ptr->vert[(n+1)%3]->x - op_ptr->x)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) 
			+ (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y));

		B = (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba - Bb);
		By= (ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(Ba_y - Bb_y);

		Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2);
		Ca_y = -2.0*(ele_ptr->vert[(n+1)%3]->y - op_ptr->y);
		Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2);
		Cb_y = -2.0*(ele_ptr->vert[n]->y - op_ptr->y);

		C = 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(log(Ca) - log(Cb));
		Cy= 0.5*(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ca_y/Ca - Cb_y/Cb);

		sum += (Ay*(B+C) + A*(By+Cy));
	}

	// 按照重力正演的传统 这里取负y方向的梯度
	return 2.0e+8*GCTL_G0*sum;
}