266 lines
9.4 KiB
C++
266 lines
9.4 KiB
C++
|
/********************************************************
|
|||
|
|
* ██████╗ ██████╗████████╗██╗
|
|||
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|||
|
|
* ██║ ███╗██║ ██║ ██║
|
|||
|
|
* ██║ ██║██║ ██║ ██║
|
|||
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|||
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|||
|
|
* 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_polygon.h"
|
|||
|
|
#include "cmath"
|
|||
|
|
|
|||
|
|
|
|||
|
|
typedef void (*gobser_polygon_ptr)(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho);
|
|||
|
|
|
|||
|
|
void gobser_polygon_vz(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho);
|
|||
|
|
void gobser_polygon_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho);
|
|||
|
|
void gobser_polygon_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho);
|
|||
|
|
|
|||
|
|
|
|||
|
|
void gctl::gobser(array<double> &out_obs, const array<vertex2dc> &cor_vert, const array<point2dc> &obsp,
|
|||
|
|
const double &rho, gravitational_field_type_e comp_id)
|
|||
|
|
{
|
|||
|
|
gobser_polygon_ptr polygon_obser;
|
|||
|
|
switch (comp_id)
|
|||
|
|
{
|
|||
|
|
case Vz:
|
|||
|
|
polygon_obser = gobser_polygon_vz;
|
|||
|
|
break;
|
|||
|
|
case Tzx:
|
|||
|
|
polygon_obser = gobser_polygon_vzx;
|
|||
|
|
break;
|
|||
|
|
case Tzz:
|
|||
|
|
polygon_obser = gobser_polygon_vzz;
|
|||
|
|
break;
|
|||
|
|
default:
|
|||
|
|
polygon_obser = gobser_polygon_vz;
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
return polygon_obser(out_obs, cor_vert, obsp, rho);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double gkernel_polygon_vz_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const gctl::point2dc &op_ptr);
|
|||
|
|
double gkernel_polygon_vzx_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const gctl::point2dc &op_ptr);
|
|||
|
|
double gkernel_polygon_vzz_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const 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 gobser_polygon_vz(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho)
|
|||
|
|
{
|
|||
|
|
int o_size = obsp.size();
|
|||
|
|
out_obs.resize(o_size, 0.0);
|
|||
|
|
|
|||
|
|
int i;
|
|||
|
|
#pragma omp parallel for private (i) shared (cor_vert, rho) schedule(guided)
|
|||
|
|
for (i = 0; i < o_size; i++)
|
|||
|
|
{
|
|||
|
|
out_obs[i] = gkernel_polygon_vz_sig(cor_vert, obsp[i]) * rho;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void gobser_polygon_vzx(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho)
|
|||
|
|
{
|
|||
|
|
int o_size = obsp.size();
|
|||
|
|
out_obs.resize(o_size, 0.0);
|
|||
|
|
|
|||
|
|
int i;
|
|||
|
|
#pragma omp parallel for private (i) shared (cor_vert, rho) schedule(guided)
|
|||
|
|
for (i = 0; i < o_size; i++)
|
|||
|
|
{
|
|||
|
|
out_obs[i] = gkernel_polygon_vzx_sig(cor_vert, obsp[i]) * rho;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void gobser_polygon_vzz(gctl::array<double> &out_obs, const gctl::array<gctl::vertex2dc> &cor_vert,
|
|||
|
|
const gctl::array<gctl::point2dc> &obsp, const double &rho)
|
|||
|
|
{
|
|||
|
|
int o_size = obsp.size();
|
|||
|
|
out_obs.resize(o_size, 0.0);
|
|||
|
|
|
|||
|
|
int i;
|
|||
|
|
#pragma omp parallel for private (i) shared (cor_vert, rho) schedule(guided)
|
|||
|
|
for (i = 0; i < o_size; i++)
|
|||
|
|
{
|
|||
|
|
out_obs[i] = gkernel_polygon_vzz_sig(cor_vert, obsp[i]) * rho;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double gkernel_polygon_vz_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const gctl::point2dc &op_ptr)
|
|||
|
|
{
|
|||
|
|
double sum;
|
|||
|
|
double A, B, C;
|
|||
|
|
double Aa, Ab, Ba, Bb, Ca, Cb;
|
|||
|
|
|
|||
|
|
int n1;
|
|||
|
|
int cor_num = cor_vert.size();
|
|||
|
|
|
|||
|
|
sum = 0.0;
|
|||
|
|
for (int n = 0; n < cor_num; n++)
|
|||
|
|
{
|
|||
|
|
n1 = (n+1)%cor_num;
|
|||
|
|
|
|||
|
|
Aa = (cor_vert[n].x - op_ptr.x)*(cor_vert[n1].y - op_ptr.y) -
|
|||
|
|
(cor_vert[n1].x - op_ptr.x)*(cor_vert[n].y - op_ptr.y);
|
|||
|
|
Ab = pow(cor_vert[n1].y - cor_vert[n].y,2) +
|
|||
|
|
pow(cor_vert[n1].x - cor_vert[n].x,2);
|
|||
|
|
A = Aa/Ab;
|
|||
|
|
|
|||
|
|
Ba = atan2(cor_vert[n].y - op_ptr.y, cor_vert[n].x - op_ptr.x);
|
|||
|
|
Bb = atan2(cor_vert[n1].y - op_ptr.y, cor_vert[n1].x - op_ptr.x);
|
|||
|
|
B = (cor_vert[n1].x - cor_vert[n].x)*(Ba - Bb);
|
|||
|
|
|
|||
|
|
Ca = pow(cor_vert[n1].x - op_ptr.x,2) + pow(cor_vert[n1].y - op_ptr.y,2);
|
|||
|
|
Cb = pow(cor_vert[n].x - op_ptr.x,2) + pow(cor_vert[n].y - op_ptr.y,2);
|
|||
|
|
C = 0.5*(cor_vert[n1].y - cor_vert[n].y)*(log(Ca) - log(Cb));
|
|||
|
|
|
|||
|
|
sum += A*(B+C);
|
|||
|
|
}
|
|||
|
|
return -2.0e+8*GCTL_G0*sum;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double gkernel_polygon_vzx_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const 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;
|
|||
|
|
|
|||
|
|
int n1;
|
|||
|
|
int cor_num = cor_vert.size();
|
|||
|
|
|
|||
|
|
sum = 0.0;
|
|||
|
|
for (int n = 0; n < cor_num; n++)
|
|||
|
|
{
|
|||
|
|
n1 = (n+1)%cor_num;
|
|||
|
|
|
|||
|
|
Aa = (cor_vert[n].x - op_ptr.x)*(cor_vert[n1].y - op_ptr.y) -
|
|||
|
|
(cor_vert[n1].x - op_ptr.x)*(cor_vert[n].y - op_ptr.y);
|
|||
|
|
Aa_x = cor_vert[n].y - cor_vert[n1].y;
|
|||
|
|
Ab = pow(cor_vert[n1].y - cor_vert[n].y,2) +
|
|||
|
|
pow(cor_vert[n1].x - cor_vert[n].x,2);
|
|||
|
|
|
|||
|
|
A = Aa/Ab;
|
|||
|
|
Ax= Aa_x/Ab;
|
|||
|
|
|
|||
|
|
Ba = atan2(cor_vert[n].y - op_ptr.y, cor_vert[n].x - op_ptr.x);
|
|||
|
|
Ba_x = (cor_vert[n].y - op_ptr.y)/((cor_vert[n].x - op_ptr.x)*(cor_vert[n].x - op_ptr.x)
|
|||
|
|
+ (cor_vert[n].y - op_ptr.y)*(cor_vert[n].y - op_ptr.y));
|
|||
|
|
|
|||
|
|
Bb = atan2(cor_vert[n1].y - op_ptr.y, cor_vert[n1].x - op_ptr.x);
|
|||
|
|
Bb_x = (cor_vert[n1].y - op_ptr.y)/((cor_vert[n1].x - op_ptr.x)*(cor_vert[n1].x - op_ptr.x)
|
|||
|
|
+ (cor_vert[n1].y - op_ptr.y)*(cor_vert[n1].y - op_ptr.y));
|
|||
|
|
|
|||
|
|
B = (cor_vert[n1].x - cor_vert[n].x)*(Ba - Bb);
|
|||
|
|
Bx= (cor_vert[n1].x - cor_vert[n].x)*(Ba_x - Bb_x);
|
|||
|
|
|
|||
|
|
Ca = pow(cor_vert[n1].x - op_ptr.x,2) + pow(cor_vert[n1].y - op_ptr.y,2);
|
|||
|
|
Ca_x = -2.0*(cor_vert[n1].x - op_ptr.x);
|
|||
|
|
Cb = pow(cor_vert[n].x - op_ptr.x,2) + pow(cor_vert[n].y - op_ptr.y,2);
|
|||
|
|
Cb_x = -2.0*(cor_vert[n].x - op_ptr.x);
|
|||
|
|
|
|||
|
|
C = 0.5*(cor_vert[n1].y - cor_vert[n].y)*(log(Ca) - log(Cb));
|
|||
|
|
Cx= 0.5*(cor_vert[n1].y - cor_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_polygon_vzz_sig(const gctl::array<gctl::vertex2dc> &cor_vert, const 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;
|
|||
|
|
|
|||
|
|
int n1;
|
|||
|
|
int cor_num = cor_vert.size();
|
|||
|
|
|
|||
|
|
sum = 0.0;
|
|||
|
|
for (int n = 0; n < cor_num; n++)
|
|||
|
|
{
|
|||
|
|
n1 = (n+1)%cor_num;
|
|||
|
|
|
|||
|
|
Aa = (cor_vert[n].x - op_ptr.x)*(cor_vert[n1].y - op_ptr.y) -
|
|||
|
|
(cor_vert[n1].x - op_ptr.x)*(cor_vert[n].y - op_ptr.y);
|
|||
|
|
Aa_y = cor_vert[n1].x - cor_vert[n].x;
|
|||
|
|
Ab = pow(cor_vert[n1].y - cor_vert[n].y,2) +
|
|||
|
|
pow(cor_vert[n1].x - cor_vert[n].x,2);
|
|||
|
|
|
|||
|
|
A = Aa/Ab;
|
|||
|
|
Ay= Aa_y/Ab;
|
|||
|
|
|
|||
|
|
Ba = atan2(cor_vert[n].y - op_ptr.y, cor_vert[n].x - op_ptr.x);
|
|||
|
|
Ba_y = -1.0*(cor_vert[n].x - op_ptr.x)/((cor_vert[n].x - op_ptr.x)*(cor_vert[n].x - op_ptr.x)
|
|||
|
|
+ (cor_vert[n].y - op_ptr.y)*(cor_vert[n].y - op_ptr.y));
|
|||
|
|
|
|||
|
|
Bb = atan2(cor_vert[n1].y - op_ptr.y, cor_vert[n1].x - op_ptr.x);
|
|||
|
|
Bb_y = -1.0*(cor_vert[n1].x - op_ptr.x)/((cor_vert[n1].x - op_ptr.x)*(cor_vert[n1].x - op_ptr.x)
|
|||
|
|
+ (cor_vert[n1].y - op_ptr.y)*(cor_vert[n1].y - op_ptr.y));
|
|||
|
|
|
|||
|
|
B = (cor_vert[n1].x - cor_vert[n].x)*(Ba - Bb);
|
|||
|
|
By= (cor_vert[n1].x - cor_vert[n].x)*(Ba_y - Bb_y);
|
|||
|
|
|
|||
|
|
Ca = pow(cor_vert[n1].x - op_ptr.x,2) + pow(cor_vert[n1].y - op_ptr.y,2);
|
|||
|
|
Ca_y = -2.0*(cor_vert[n1].y - op_ptr.y);
|
|||
|
|
Cb = pow(cor_vert[n].x - op_ptr.x,2) + pow(cor_vert[n].y - op_ptr.y,2);
|
|||
|
|
Cb_y = -2.0*(cor_vert[n].y - op_ptr.y);
|
|||
|
|
|
|||
|
|
C = 0.5*(cor_vert[n1].y - cor_vert[n].y)*(log(Ca) - log(Cb));
|
|||
|
|
Cy= 0.5*(cor_vert[n1].y - cor_vert[n].y)*(Ca_y/Ca - Cb_y/Cb);
|
|||
|
|
|
|||
|
|
sum += (Ay*(B+C) + A*(By+Cy));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// 按照重力正演的传统 这里取负y方向的梯度
|
|||
|
|
return 2.0e+8*GCTL_G0*sum;
|
|||
|
|
}
|