/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * 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 . * * 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 &out_obs, const gctl::array &cor_vert, const gctl::array &obsp, const double &rho); void gobser_polygon_vz(gctl::array &out_obs, const gctl::array &cor_vert, const gctl::array &obsp, const double &rho); void gobser_polygon_vzx(gctl::array &out_obs, const gctl::array &cor_vert, const gctl::array &obsp, const double &rho); void gobser_polygon_vzz(gctl::array &out_obs, const gctl::array &cor_vert, const gctl::array &obsp, const double &rho); void gctl::gobser(array &out_obs, const array &cor_vert, const array &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 &cor_vert, const gctl::point2dc &op_ptr); double gkernel_polygon_vzx_sig(const gctl::array &cor_vert, const gctl::point2dc &op_ptr); double gkernel_polygon_vzz_sig(const gctl::array &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 &out_obs, const gctl::array &cor_vert, const gctl::array &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 &out_obs, const gctl::array &cor_vert, const gctl::array &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 &out_obs, const gctl::array &cor_vert, const gctl::array &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 &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 &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 &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; }