zhangyiss/gctl_potential
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
gctl_potential/lib/potential/mkernel_block.cpp
Yi Zhang 042d310021 initial upload
2024-09-10 19:56:41 +08:00

1122 lines
53 KiB
C++
Raw Permalink Blame History

/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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_block.h"
#include "cmath"
/**
* @brief callback interface of the magnetic kernel of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param geo_inclina Inclination angle of the geo-magnetic field.
* @param geo_declina Declination angle of the geo-magnetic field.
* @param[in] comp_id The component of magnetic field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an matrix<double> object
*/
typedef void (*magkernel_block_ptr)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, gctl::verbose_type_e verbose);
void magkernel_block_deltaT(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, gctl::verbose_type_e verbose);
void magkernel_block_deltaTx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, gctl::verbose_type_e verbose);
void magkernel_block_deltaTy(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, gctl::verbose_type_e verbose);
void magkernel_block_deltaTz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, gctl::verbose_type_e verbose);
/**
* @brief Integrated callback function of the magnetic kernel of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param geo_inclina Inclination angle of the geo-magnetic field.
* @param geo_declinae Declination angle of the geo-magnetic field.
* @param[in] comp_id The component of magnetic field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an matrix<double> object
*/
void gctl::magkernel(matrix<double> &out_kernel, const array<mag_block> &ele, const array<point3dc> &obsp,
double geo_inclina, double geo_declina, magnetic_field_type_e comp_id, verbose_type_e verbose)
{
magkernel_block_ptr block_kernel;
switch (comp_id)
{
case DeltaT:
block_kernel = magkernel_block_deltaT;
break;
case DeltaTx:
block_kernel = magkernel_block_deltaTx;
break;
case DeltaTy:
block_kernel = magkernel_block_deltaTy;
break;
case DeltaTz:
block_kernel = magkernel_block_deltaTz;
break;
default:
block_kernel = magkernel_block_deltaT;
break;
}
return block_kernel(out_kernel, ele, obsp, geo_inclina, geo_declina, verbose);
}
/**
* @brief callback interface of the magnetic kernel of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] comp_id The component of magnetic field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an matrix<double> object
*/
typedef void (*magkernel_block_ptr2)(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void magkernel_block_hax(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void magkernel_block_hay(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
void magkernel_block_za(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, gctl::verbose_type_e verbose);
/**
* @brief Integrated callback function of the magnetic kernel of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] comp_id The component of magnetic field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an matrix<double> object
*/
void gctl::magkernel(matrix<double> &out_kernel, const array<mag_block> &ele,
const array<point3dc> &obsp, magnetic_field_type_e comp_id, verbose_type_e verbose)
{
magkernel_block_ptr2 block_kernel;
switch (comp_id)
{
case Hax:
block_kernel = magkernel_block_hax;
break;
case Hay:
block_kernel = magkernel_block_hay;
break;
case Za:
block_kernel = magkernel_block_za;
break;
default:
block_kernel = magkernel_block_za;
break;
}
return block_kernel(out_kernel, ele, obsp, verbose);
}
/**
* @brief callback interface of the gravitational observation of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] mag_ptr Magnetism of the calculating elements
* @param geo_inclina Inclination angle of the geo-magnetic field.
* @param geo_declina Declination angle of the geo-magnetic field.
* @param[in] comp_id The component of gravitational field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an array<double> object
*/
typedef void (*magobser_block_ptr)(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, gctl::verbose_type_e verbose);
void magobser_block_deltaT(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, gctl::verbose_type_e verbose);
void magobser_block_deltaTx(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, gctl::verbose_type_e verbose);
void magobser_block_deltaTy(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, gctl::verbose_type_e verbose);
void magobser_block_deltaTz(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, gctl::verbose_type_e verbose);
/**
* @brief Integrated callback function of the gravitational observation of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] mag_ptr Magnetism of the calculating elements
* @param geo_inclina Inclination angle of the geo-magnetic field.
* @param geo_declina Declination angle of the geo-magnetic field.
* @param[in] comp_id The component of gravitational field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an array<double> object
*/
void gctl::magobser(array<double> &out_obs, const array<mag_block> &ele, const array<point3dc> &obsp,
const array<double> &sus, double geo_inclina, double geo_declina, magnetic_field_type_e comp_id,
verbose_type_e verbose)
{
magobser_block_ptr block_obser;
switch (comp_id)
{
case DeltaT:
block_obser = magobser_block_deltaT;
break;
case DeltaTx:
block_obser = magobser_block_deltaTx;
break;
case DeltaTy:
block_obser = magobser_block_deltaTy;
break;
case DeltaTz:
block_obser = magobser_block_deltaTz;
break;
default:
block_obser = magobser_block_deltaT;
break;
}
return block_obser(out_obs, ele, obsp, sus, geo_inclina, geo_declina, verbose);
}
/**
* @brief callback interface of the gravitational observation of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] mag_ptr Magnetism of the calculating elements
* @param[in] comp_id The component of gravitational field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an array<double> object
*/
typedef void (*magobser_block_ptr2)(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, gctl::verbose_type_e verbose);
void magobser_block_hax(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, gctl::verbose_type_e verbose);
void magobser_block_hay(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, gctl::verbose_type_e verbose);
void magobser_block_za(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, gctl::verbose_type_e verbose);
/**
* @brief Integrated callback function of the gravitational observation of given elements type.
*
* @param[in] ele_ptr The element's pointer.
* @param op_ptr Coordinates of the observation.
* @param[in] mag_ptr Magnetism of the calculating elements
* @param geo_inclina Inclination angle of the geo-magnetic field.
* @param geo_declina Declination angle of the geo-magnetic field.
* @param[in] comp_id The component of gravitational field needs to be calculated.
* @param[in] show_progress The show progress
*
* @return Pointer of an array<double> object
*/
void gctl::magobser(array<double> &out_obs, const array<mag_block> &ele, const array<point3dc> &obsp,
const array<double> &sus, magnetic_field_type_e comp_id, verbose_type_e verbose)
{
magobser_block_ptr2 block_obser;
switch (comp_id)
{
case Hax:
block_obser = magobser_block_hax;
break;
case Hay:
block_obser = magobser_block_hay;
break;
case Za:
block_obser = magobser_block_za;
break;
default:
block_obser = magobser_block_za;
break;
}
return block_obser(out_obs, ele, obsp, sus, verbose);
}
// declare algorithm for individual element and observation point
double magkernel_block_deltaT_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina, double geo_declina);
double magkernel_block_deltaTx_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina, double geo_declina);
double magkernel_block_deltaTy_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina, double geo_declina);
double magkernel_block_deltaTz_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina, double geo_declina);
double magkernel_block_hax_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr);
double magkernel_block_hay_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr);
double magkernel_block_za_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr);
// define functions
void magkernel_block_deltaT(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, 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, "magkernel_deltaT");
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] = magkernel_block_deltaT_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina);
}
}
return;
}
void magkernel_block_deltaTx(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, 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, "magkernel_deltaTx");
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] = magkernel_block_deltaTx_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina);
}
}
return;
}
void magkernel_block_deltaTy(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, 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, "magkernel_deltaTy");
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] = magkernel_block_deltaTy_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina);
}
}
return;
}
void magkernel_block_deltaTz(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, double geo_inclina, double geo_declina, 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, "magkernel_deltaTz");
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] = magkernel_block_deltaTz_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina);
}
}
return;
}
void magkernel_block_hax(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &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, "magkernel_Hax");
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] = magkernel_block_hax_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void magkernel_block_hay(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &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, "magkernel_Hay");
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] = magkernel_block_hay_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void magkernel_block_za(gctl::matrix<double> &out_kernel, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &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, "magkernel_Za");
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] = magkernel_block_za_sig(ele.get(j), obsp.get(i));
}
}
return;
}
void magobser_block_deltaT(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, 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, "magobser_deltaT");
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] += magkernel_block_deltaT_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina) * sus[j];
}
}
return;
}
void magobser_block_deltaTx(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, 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, "magobser_deltaTx");
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] += magkernel_block_deltaTx_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina) * sus[j];
}
}
return;
}
void magobser_block_deltaTy(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, 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, "magobser_deltaTy");
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] += magkernel_block_deltaTy_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina) * sus[j];
}
}
return;
}
void magobser_block_deltaTz(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, double geo_inclina,
double geo_declina, 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, "magobser_deltaTz");
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] += magkernel_block_deltaTz_sig(ele.get(j), obsp.get(i), geo_inclina, geo_declina) * sus[j];
}
}
return;
}
void magobser_block_hax(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, 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, "magobser_Hax");
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] += magkernel_block_hax_sig(ele.get(j), obsp.get(i)) * sus[j];
}
}
return;
}
void magobser_block_hay(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, 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, "magobser_Hay");
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] += magkernel_block_hay_sig(ele.get(j), obsp.get(i)) * sus[j];
}
}
return;
}
void magobser_block_za(gctl::array<double> &out_obs, const gctl::array<gctl::mag_block> &ele,
const gctl::array<gctl::point3dc> &obsp, const gctl::array<double> &sus, 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, "magobser_Za");
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] += magkernel_block_za_sig(ele.get(j), obsp.get(i)) * sus[j];
}
}
return;
}
double magkernel_block_deltaT_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina,
double geo_declina)
{
double I0,A0,I,A;
double k1,k2,k3,k4,k5,k6;
double x1,x2,y1,y2,z1,z2;
double w111,w112,w121,w122,w211,w212,w221,w222;
double v111,v112,v121,v122,v211,v212,v221,v222;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I0= geo_inclina*GCTL_Pi/180.0;
A0= (90 - geo_declina)*GCTL_Pi/180.0;
I = (mp->inclina_deg)*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
k4=cos(I0)*cos(A0)*cos(I)*cos(A);
k5=cos(I0)*sin(A0)*cos(I)*sin(A);
k6=-sin(I0)*sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222=((x2-op_ptr->x)*(x2-op_ptr->x))+R222*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
w122=((x1-op_ptr->x)*(x1-op_ptr->x))+R122*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
w212=((x2-op_ptr->x)*(x2-op_ptr->x))+R212*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
w112=((x1-op_ptr->x)*(x1-op_ptr->x))+R112*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
w221=((x2-op_ptr->x)*(x2-op_ptr->x))+R221*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
w121=((x1-op_ptr->x)*(x1-op_ptr->x))+R121*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
w211=((x2-op_ptr->x)*(x2-op_ptr->x))+R211*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
w111=((x1-op_ptr->x)*(x1-op_ptr->x))+R111*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
v222=((y2-op_ptr->y)*(y2-op_ptr->y))+R222*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
v122=((y2-op_ptr->y)*(y2-op_ptr->y))+R122*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
v212=((y1-op_ptr->y)*(y1-op_ptr->y))+R212*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
v112=((y1-op_ptr->y)*(y1-op_ptr->y))+R112*(z2-op_ptr->z)+((z2-op_ptr->z)*(z2-op_ptr->z));
v221=((y2-op_ptr->y)*(y2-op_ptr->y))+R221*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
v121=((y2-op_ptr->y)*(y2-op_ptr->y))+R121*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
v211=((y1-op_ptr->y)*(y1-op_ptr->y))+R211*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
v111=((y1-op_ptr->y)*(y1-op_ptr->y))+R111*(z1-op_ptr->z)+((z1-op_ptr->z)*(z1-op_ptr->z));
G222=k1*log(R222+x2-op_ptr->x)+k2*log(R222+y2-op_ptr->y)+k3*log(R222+z2-op_ptr->z)
+k4*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/w222)+k5*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/v222)
+k6*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/(R222*(z2-op_ptr->z)));
G122=k1*log(R122+x1-op_ptr->x)+k2*log(R122+y2-op_ptr->y)+k3*log(R122+z2-op_ptr->z)
+k4*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/w122)+k5*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/v122)
+k6*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/(R122*(z2-op_ptr->z)));
G212=k1*log(R212+x2-op_ptr->x)+k2*log(R212+y1-op_ptr->y)+k3*log(R212+z2-op_ptr->z)
+k4*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/w212)+k5*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/v212)
+k6*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/(R212*(z2-op_ptr->z)));
G112=k1*log(R112+x1-op_ptr->x)+k2*log(R112+y1-op_ptr->y)+k3*log(R112+z2-op_ptr->z)
+k4*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/w112)+k5*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/v112)
+k6*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/(R112*(z2-op_ptr->z)));
G221=k1*log(R221+x2-op_ptr->x)+k2*log(R221+y2-op_ptr->y)+k3*log(R221+z1-op_ptr->z)
+k4*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/w221)+k5*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/v221)
+k6*gctl::arctg(((x2-op_ptr->x)*(y2-op_ptr->y))/(R221*(z1-op_ptr->z)));
G121=k1*log(R121+x1-op_ptr->x)+k2*log(R121+y2-op_ptr->y)+k3*log(R121+z1-op_ptr->z)
+k4*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/w121)+k5*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/v121)
+k6*gctl::arctg(((x1-op_ptr->x)*(y2-op_ptr->y))/(R121*(z1-op_ptr->z)));
G211=k1*log(R211+x2-op_ptr->x)+k2*log(R211+y1-op_ptr->y)+k3*log(R211+z1-op_ptr->z)
+k4*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/w211)+k5*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/v211)
+k6*gctl::arctg(((x2-op_ptr->x)*(y1-op_ptr->y))/(R211*(z1-op_ptr->z)));
G111=k1*log(R111+x1-op_ptr->x)+k2*log(R111+y1-op_ptr->y)+k3*log(R111+z1-op_ptr->z)
+k4*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/w111)+k5*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/v111)
+k6*gctl::arctg(((x1-op_ptr->x)*(y1-op_ptr->y))/(R111*(z1-op_ptr->z)));
return GCTL_T0/(4.0*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_deltaTx_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina,
double geo_declina)
{
double I0,A0,I,A;
double k1,k2,k3,k4,k5,k6;
double x1,x2,y1,y2,z1,z2;
double w111,w112,w121,w122,w211,w212,w221,w222;
double v111,v112,v121,v122,v211,v212,v221,v222;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I0= geo_inclina*GCTL_Pi/180.0;
A0= (90 - geo_declina)*GCTL_Pi/180.0;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
k4=cos(I0)*cos(A0)*cos(I)*cos(A);
k5=cos(I0)*sin(A0)*cos(I)*sin(A);
k6=-sin(I0)*sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222=((y2-op_ptr->y)*(pow(x2-op_ptr->x,2)-R222*(z2-op_ptr->z)))/(R222*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2))*(R222+z2-op_ptr->z));
w122=((y2-op_ptr->y)*(pow(x1-op_ptr->x,2)-R122*(z2-op_ptr->z)))/(R122*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2))*(R122+z2-op_ptr->z));
w212=((y1-op_ptr->y)*(pow(x2-op_ptr->x,2)-R212*(z2-op_ptr->z)))/(R212*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2))*(R212+z2-op_ptr->z));
w112=((y1-op_ptr->y)*(pow(x1-op_ptr->x,2)-R112*(z2-op_ptr->z)))/(R112*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2))*(R112+z2-op_ptr->z));
w221=((y2-op_ptr->y)*(pow(x2-op_ptr->x,2)-R221*(z1-op_ptr->z)))/(R221*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2))*(R221+z1-op_ptr->z));
w121=((y2-op_ptr->y)*(pow(x1-op_ptr->x,2)-R121*(z1-op_ptr->z)))/(R121*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2))*(R121+z1-op_ptr->z));
w211=((y1-op_ptr->y)*(pow(x2-op_ptr->x,2)-R211*(z1-op_ptr->z)))/(R211*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2))*(R211+z1-op_ptr->z));
w111=((y1-op_ptr->y)*(pow(x1-op_ptr->x,2)-R111*(z1-op_ptr->z)))/(R111*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2))*(R111+z1-op_ptr->z));
v222=((y2-op_ptr->y)*(z2-op_ptr->z))/(R222*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
v122=((y2-op_ptr->y)*(z2-op_ptr->z))/(R122*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
v212=((y1-op_ptr->y)*(z2-op_ptr->z))/(R212*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
v112=((y1-op_ptr->y)*(z2-op_ptr->z))/(R112*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
v221=((y2-op_ptr->y)*(z1-op_ptr->z))/(R221*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
v121=((y2-op_ptr->y)*(z1-op_ptr->z))/(R121*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
v211=((y1-op_ptr->y)*(z1-op_ptr->z))/(R211*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
v111=((y1-op_ptr->y)*(z1-op_ptr->z))/(R111*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
G222=-k1*pow(R222,-1)-k2*(x2-op_ptr->x)/(R222*(R222+y2-op_ptr->y))-k3*(x2-op_ptr->x)/(R222*(R222+z2-op_ptr->z))+k4*w222-k5*(y2-op_ptr->y)/(R222*(R222+z2-op_ptr->z))-k6*v222;
G122=-k1*pow(R122,-1)-k2*(x1-op_ptr->x)/(R122*(R122+y2-op_ptr->y))-k3*(x1-op_ptr->x)/(R122*(R122+z2-op_ptr->z))+k4*w122-k5*(y2-op_ptr->y)/(R122*(R122+z2-op_ptr->z))-k6*v122;
G212=-k1*pow(R212,-1)-k2*(x2-op_ptr->x)/(R212*(R212+y1-op_ptr->y))-k3*(x2-op_ptr->x)/(R212*(R212+z2-op_ptr->z))+k4*w212-k5*(y1-op_ptr->y)/(R212*(R212+z2-op_ptr->z))-k6*v212;
G112=-k1*pow(R112,-1)-k2*(x1-op_ptr->x)/(R112*(R112+y1-op_ptr->y))-k3*(x1-op_ptr->x)/(R112*(R112+z2-op_ptr->z))+k4*w112-k5*(y1-op_ptr->y)/(R112*(R112+z2-op_ptr->z))-k6*v112;
G221=-k1*pow(R221,-1)-k2*(x2-op_ptr->x)/(R221*(R221+y2-op_ptr->y))-k3*(x2-op_ptr->x)/(R221*(R221+z1-op_ptr->z))+k4*w221-k5*(y2-op_ptr->y)/(R221*(R221+z1-op_ptr->z))-k6*v221;
G121=-k1*pow(R121,-1)-k2*(x1-op_ptr->x)/(R121*(R121+y2-op_ptr->y))-k3*(x1-op_ptr->x)/(R121*(R121+z1-op_ptr->z))+k4*w121-k5*(y2-op_ptr->y)/(R121*(R121+z1-op_ptr->z))-k6*v121;
G211=-k1*pow(R211,-1)-k2*(x2-op_ptr->x)/(R211*(R211+y1-op_ptr->y))-k3*(x2-op_ptr->x)/(R211*(R211+z1-op_ptr->z))+k4*w211-k5*(y1-op_ptr->y)/(R211*(R211+z1-op_ptr->z))-k6*v211;
G111=-k1*pow(R111,-1)-k2*(x1-op_ptr->x)/(R111*(R111+y1-op_ptr->y))-k3*(x1-op_ptr->x)/(R111*(R111+z1-op_ptr->z))+k4*w111-k5*(y1-op_ptr->y)/(R111*(R111+z1-op_ptr->z))-k6*v111;
return GCTL_T0/(4.0*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_deltaTy_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina,
double geo_declina)
{
double I0,A0,I,A;
double k1,k2,k3,k4,k5,k6;
double x1,x2,y1,y2,z1,z2;
double w111,w112,w121,w122,w211,w212,w221,w222;
double v111,v112,v121,v122,v211,v212,v221,v222;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I0= geo_inclina*GCTL_Pi/180.0;
A0= (90 - geo_declina)*GCTL_Pi/180.0;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
k4=cos(I0)*cos(A0)*cos(I)*cos(A);
k5=cos(I0)*sin(A0)*cos(I)*sin(A);
k6=-sin(I0)*sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222=((x2-op_ptr->x)*(pow(y2-op_ptr->y,2)-R222*(z2-op_ptr->z)))/(R222*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(R222+z2-op_ptr->z));
w122=((x1-op_ptr->x)*(pow(y2-op_ptr->y,2)-R122*(z2-op_ptr->z)))/(R122*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(R122+z2-op_ptr->z));
w212=((x2-op_ptr->x)*(pow(y1-op_ptr->y,2)-R212*(z2-op_ptr->z)))/(R212*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(R212+z2-op_ptr->z));
w112=((x1-op_ptr->x)*(pow(y1-op_ptr->y,2)-R112*(z2-op_ptr->z)))/(R112*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(R112+z2-op_ptr->z));
w221=((x2-op_ptr->x)*(pow(y2-op_ptr->y,2)-R221*(z1-op_ptr->z)))/(R221*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(R221+z1-op_ptr->z));
w121=((x1-op_ptr->x)*(pow(y2-op_ptr->y,2)-R121*(z1-op_ptr->z)))/(R121*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(R121+z1-op_ptr->z));
w211=((x2-op_ptr->x)*(pow(y1-op_ptr->y,2)-R211*(z1-op_ptr->z)))/(R211*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(R211+z1-op_ptr->z));
w111=((x1-op_ptr->x)*(pow(y1-op_ptr->y,2)-R111*(z1-op_ptr->z)))/(R111*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(R111+z1-op_ptr->z));
v222=((x2-op_ptr->x)*(z2-op_ptr->z))/(R222*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v122=((x1-op_ptr->x)*(z2-op_ptr->z))/(R122*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v212=((x2-op_ptr->x)*(z2-op_ptr->z))/(R212*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v112=((x1-op_ptr->x)*(z2-op_ptr->z))/(R112*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v221=((x2-op_ptr->x)*(z1-op_ptr->z))/(R221*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v121=((x1-op_ptr->x)*(z1-op_ptr->z))/(R121*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v211=((x2-op_ptr->x)*(z1-op_ptr->z))/(R211*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v111=((x1-op_ptr->x)*(z1-op_ptr->z))/(R111*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
G222=-k1*(y2-op_ptr->y)/(R222*(R222+x2-op_ptr->x))-k2*pow(R222,-1)-k3*(y2-op_ptr->y)/(R222*(R222+z2-op_ptr->z))-k4*(x2-op_ptr->x)/(R222*(R222+z2-op_ptr->z))+k5*w222-k6*v222;
G122=-k1*(y2-op_ptr->y)/(R122*(R122+x1-op_ptr->x))-k2*pow(R122,-1)-k3*(y2-op_ptr->y)/(R122*(R122+z2-op_ptr->z))-k4*(x1-op_ptr->x)/(R122*(R122+z2-op_ptr->z))+k5*w122-k6*v122;
G212=-k1*(y1-op_ptr->y)/(R212*(R212+x2-op_ptr->x))-k2*pow(R212,-1)-k3*(y1-op_ptr->y)/(R212*(R212+z2-op_ptr->z))-k4*(x2-op_ptr->x)/(R212*(R212+z2-op_ptr->z))+k5*w212-k6*v212;
G112=-k1*(y1-op_ptr->y)/(R112*(R112+x1-op_ptr->x))-k2*pow(R112,-1)-k3*(y1-op_ptr->y)/(R112*(R112+z2-op_ptr->z))-k4*(x1-op_ptr->x)/(R112*(R112+z2-op_ptr->z))+k5*w112-k6*v112;
G221=-k1*(y2-op_ptr->y)/(R221*(R221+x2-op_ptr->x))-k2*pow(R221,-1)-k3*(y2-op_ptr->y)/(R221*(R221+z1-op_ptr->z))-k4*(x2-op_ptr->x)/(R221*(R221+z1-op_ptr->z))+k5*w221-k6*v221;
G121=-k1*(y2-op_ptr->y)/(R121*(R121+x1-op_ptr->x))-k2*pow(R121,-1)-k3*(y2-op_ptr->y)/(R121*(R121+z1-op_ptr->z))-k4*(x1-op_ptr->x)/(R121*(R121+z1-op_ptr->z))+k5*w121-k6*v121;
G211=-k1*(y1-op_ptr->y)/(R211*(R211+x2-op_ptr->x))-k2*pow(R211,-1)-k3*(y1-op_ptr->y)/(R211*(R211+z1-op_ptr->z))-k4*(x2-op_ptr->x)/(R211*(R211+z1-op_ptr->z))+k5*w211-k6*v211;
G111=-k1*(y1-op_ptr->y)/(R111*(R111+x1-op_ptr->x))-k2*pow(R111,-1)-k3*(y1-op_ptr->y)/(R111*(R111+z1-op_ptr->z))-k4*(x1-op_ptr->x)/(R111*(R111+z1-op_ptr->z))+k5*w111-k6*v111;
return GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_deltaTz_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr, double geo_inclina,
double geo_declina)
{
double I0,A0,I,A;
double k1,k2,k3,k4,k5,k6;
double x1,x2,y1,y2,z1,z2;
double w111,w112,w121,w122,w211,w212,w221,w222;
double v111,v112,v121,v122,v211,v212,v221,v222;
double t111,t112,t121,t122,t211,t212,t221,t222;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I0= geo_inclina*GCTL_Pi/180.0;
A0= (90 - geo_declina)*GCTL_Pi/180.0;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
k4=cos(I0)*cos(A0)*cos(I)*cos(A);
k5=cos(I0)*sin(A0)*cos(I)*sin(A);
k6=-sin(I0)*sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222=((x2-op_ptr->x)*(y2-op_ptr->y)*(pow(z2-op_ptr->z,2)+pow(R222,2)))/(R222*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
w122=((x1-op_ptr->x)*(y2-op_ptr->y)*(pow(z2-op_ptr->z,2)+pow(R122,2)))/(R122*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
w212=((x2-op_ptr->x)*(y1-op_ptr->y)*(pow(z2-op_ptr->z,2)+pow(R212,2)))/(R212*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
w112=((x1-op_ptr->x)*(y1-op_ptr->y)*(pow(z2-op_ptr->z,2)+pow(R112,2)))/(R112*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
w221=((x2-op_ptr->x)*(y2-op_ptr->y)*(pow(z1-op_ptr->z,2)+pow(R221,2)))/(R221*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
w121=((x1-op_ptr->x)*(y2-op_ptr->y)*(pow(z1-op_ptr->z,2)+pow(R121,2)))/(R121*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
w211=((x2-op_ptr->x)*(y1-op_ptr->y)*(pow(z1-op_ptr->z,2)+pow(R211,2)))/(R211*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
w111=((x1-op_ptr->x)*(y1-op_ptr->y)*(pow(z1-op_ptr->z,2)+pow(R111,2)))/(R111*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
v222=((x2-op_ptr->x)*(y2-op_ptr->y))/(R222*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v122=((x1-op_ptr->x)*(y2-op_ptr->y))/(R122*(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v212=((x2-op_ptr->x)*(y1-op_ptr->y))/(R212*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v112=((x1-op_ptr->x)*(y1-op_ptr->y))/(R112*(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2)));
v221=((x2-op_ptr->x)*(y2-op_ptr->y))/(R221*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v121=((x1-op_ptr->x)*(y2-op_ptr->y))/(R121*(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v211=((x2-op_ptr->x)*(y1-op_ptr->y))/(R211*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
v111=((x1-op_ptr->x)*(y1-op_ptr->y))/(R111*(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2)));
t222=((x2-op_ptr->x)*(y2-op_ptr->y))/(R222*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
t122=((x1-op_ptr->x)*(y2-op_ptr->y))/(R122*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
t212=((x2-op_ptr->x)*(y1-op_ptr->y))/(R212*(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
t112=((x1-op_ptr->x)*(y1-op_ptr->y))/(R112*(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2)));
t221=((x2-op_ptr->x)*(y2-op_ptr->y))/(R221*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
t121=((x1-op_ptr->x)*(y2-op_ptr->y))/(R121*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
t211=((x2-op_ptr->x)*(y1-op_ptr->y))/(R211*(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
t111=((x1-op_ptr->x)*(y1-op_ptr->y))/(R111*(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2)));
G222=-k1*(z2-op_ptr->z)/(R222*(R222+x2-op_ptr->x))-k2*(z2-op_ptr->z)/(R222*(R222+y2-op_ptr->y))-k3/R222+k4*t222+k5*v222+k6*w222;
G122=-k1*(z2-op_ptr->z)/(R122*(R122+x1-op_ptr->x))-k2*(z2-op_ptr->z)/(R122*(R122+y2-op_ptr->y))-k3/R122+k4*t122+k5*v122+k6*w122;
G212=-k1*(z2-op_ptr->z)/(R212*(R212+x2-op_ptr->x))-k2*(z2-op_ptr->z)/(R212*(R212+y1-op_ptr->y))-k3/R212+k4*t212+k5*v212+k6*w212;
G112=-k1*(z2-op_ptr->z)/(R112*(R112+x1-op_ptr->x))-k2*(z2-op_ptr->z)/(R112*(R112+y1-op_ptr->y))-k3/R112+k4*t112+k5*v112+k6*w112;
G221=-k1*(z1-op_ptr->z)/(R221*(R221+x2-op_ptr->x))-k2*(z1-op_ptr->z)/(R221*(R221+y2-op_ptr->y))-k3/R221+k4*t221+k5*v221+k6*w221;
G121=-k1*(z1-op_ptr->z)/(R121*(R121+x1-op_ptr->x))-k2*(z1-op_ptr->z)/(R121*(R121+y2-op_ptr->y))-k3/R121+k4*t121+k5*v121+k6*w121;
G211=-k1*(z1-op_ptr->z)/(R211*(R211+x2-op_ptr->x))-k2*(z1-op_ptr->z)/(R211*(R211+y1-op_ptr->y))-k3/R211+k4*t211+k5*v211+k6*w211;
G111=-k1*(z1-op_ptr->z)/(R111*(R111+x1-op_ptr->x))-k2*(z1-op_ptr->z)/(R111*(R111+y1-op_ptr->y))-k3/R111+k4*t111+k5*v111+k6*w111;
return GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_hax_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr)
{
double I,A;
double Alpha,Beta,Gamma;
double x1,x2,y1,y2,z1,z2;
double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
Alpha=cos(I)*cos(A);
Beta=cos(I)*sin(A);
Gamma=sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w122_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w212_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w112_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w221_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w121_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w211_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w111_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w222_2=(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2))+R222*(z2-op_ptr->z);
w122_2=(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2))+R122*(z2-op_ptr->z);
w212_2=(pow(x2-op_ptr->x,2)+pow(z2-op_ptr->z,2))+R212*(z2-op_ptr->z);
w112_2=(pow(x1-op_ptr->x,2)+pow(z2-op_ptr->z,2))+R112*(z2-op_ptr->z);
w221_2=(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2))+R221*(z1-op_ptr->z);
w121_2=(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2))+R121*(z1-op_ptr->z);
w211_2=(pow(x2-op_ptr->x,2)+pow(z1-op_ptr->z,2))+R211*(z1-op_ptr->z);
w111_2=(pow(x1-op_ptr->x,2)+pow(z1-op_ptr->z,2))+R111*(z1-op_ptr->z);
G222=Alpha*atan2(w222_1,w222_2)+Beta*log(R222+z2-op_ptr->z)+Gamma*log(R222+y2-op_ptr->y);
G122=Alpha*atan2(w122_1,w122_2)+Beta*log(R122+z2-op_ptr->z)+Gamma*log(R122+y2-op_ptr->y);
G212=Alpha*atan2(w212_1,w212_2)+Beta*log(R212+z2-op_ptr->z)+Gamma*log(R212+y1-op_ptr->y);
G112=Alpha*atan2(w112_1,w112_2)+Beta*log(R112+z2-op_ptr->z)+Gamma*log(R112+y1-op_ptr->y);
G221=Alpha*atan2(w221_1,w221_2)+Beta*log(R221+z1-op_ptr->z)+Gamma*log(R221+y2-op_ptr->y);
G121=Alpha*atan2(w121_1,w121_2)+Beta*log(R121+z1-op_ptr->z)+Gamma*log(R121+y2-op_ptr->y);
G211=Alpha*atan2(w211_1,w211_2)+Beta*log(R211+z1-op_ptr->z)+Gamma*log(R211+y1-op_ptr->y);
G111=Alpha*atan2(w111_1,w111_2)+Beta*log(R111+z1-op_ptr->z)+Gamma*log(R111+y1-op_ptr->y);
return GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_hay_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr)
{
double I,A;
double Alpha,Beta,Gamma;
double x1,x2,y1,y2,z1,z2;
double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
Alpha=cos(I)*cos(A);
Beta=cos(I)*sin(A);
Gamma=sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w122_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w212_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w112_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w221_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w121_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w211_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w111_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w222_2=(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))+R222*(z2-op_ptr->z);
w122_2=(pow(y2-op_ptr->y,2)+pow(z2-op_ptr->z,2))+R122*(z2-op_ptr->z);
w212_2=(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))+R212*(z2-op_ptr->z);
w112_2=(pow(y1-op_ptr->y,2)+pow(z2-op_ptr->z,2))+R112*(z2-op_ptr->z);
w221_2=(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))+R221*(z1-op_ptr->z);
w121_2=(pow(y2-op_ptr->y,2)+pow(z1-op_ptr->z,2))+R121*(z1-op_ptr->z);
w211_2=(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))+R211*(z1-op_ptr->z);
w111_2=(pow(y1-op_ptr->y,2)+pow(z1-op_ptr->z,2))+R111*(z1-op_ptr->z);
G222=Beta*atan2(w222_1,w222_2)+Alpha*log(R222+z2-op_ptr->z)+Gamma*log(R222+x2-op_ptr->x);
G122=Beta*atan2(w122_1,w122_2)+Alpha*log(R122+z2-op_ptr->z)+Gamma*log(R122+x1-op_ptr->x);
G212=Beta*atan2(w212_1,w212_2)+Alpha*log(R212+z2-op_ptr->z)+Gamma*log(R212+x2-op_ptr->x);
G112=Beta*atan2(w112_1,w112_2)+Alpha*log(R112+z2-op_ptr->z)+Gamma*log(R112+x1-op_ptr->x);
G221=Beta*atan2(w221_1,w221_2)+Alpha*log(R221+z1-op_ptr->z)+Gamma*log(R221+x2-op_ptr->x);
G121=Beta*atan2(w121_1,w121_2)+Alpha*log(R121+z1-op_ptr->z)+Gamma*log(R121+x1-op_ptr->x);
G211=Beta*atan2(w211_1,w211_2)+Alpha*log(R211+z1-op_ptr->z)+Gamma*log(R211+x2-op_ptr->x);
G111=Beta*atan2(w111_1,w111_2)+Alpha*log(R111+z1-op_ptr->z)+Gamma*log(R111+x1-op_ptr->x);
return GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
double magkernel_block_za_sig(gctl::mag_block *ele_ptr, gctl::point3dc *op_ptr)
{
double I,A;
double Alpha,Beta,Gamma;
double x1,x2,y1,y2,z1,z2;
double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
double R222,R122,R212,R112,R221,R121,R211,R111;
double G222,G122,G212,G112,G221,G121,G211,G111;
if (ele_ptr->att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
gctl::magblock_para* mp = ele_ptr->att;
I = mp->inclina_deg*GCTL_Pi/180.0;
A = (90 - mp->declina_deg)*GCTL_Pi/180.0;
Alpha=cos(I)*cos(A);
Beta=cos(I)*sin(A);
Gamma=sin(I);
x1 = ele_ptr->dl->x; x2 = ele_ptr->ur->x;
y1 = ele_ptr->dl->y; y2 = ele_ptr->ur->y;
z1 = ele_ptr->dl->z; z2 = ele_ptr->ur->z;
R222=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R122=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R212=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R112=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z2-op_ptr->z)*(z2-op_ptr->z));
R221=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R121=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y2-op_ptr->y)*(y2-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R211=sqrt((x2-op_ptr->x)*(x2-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
R111=sqrt((x1-op_ptr->x)*(x1-op_ptr->x)+(y1-op_ptr->y)*(y1-op_ptr->y)+(z1-op_ptr->z)*(z1-op_ptr->z));
w222_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w122_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w212_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w112_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w221_1=(x2-op_ptr->x)*(y2-op_ptr->y);
w121_1=(x1-op_ptr->x)*(y2-op_ptr->y);
w211_1=(x2-op_ptr->x)*(y1-op_ptr->y);
w111_1=(x1-op_ptr->x)*(y1-op_ptr->y);
w222_2=R222*(z2-op_ptr->z);
w122_2=R122*(z2-op_ptr->z);
w212_2=R212*(z2-op_ptr->z);
w112_2=R112*(z2-op_ptr->z);
w221_2=R221*(z1-op_ptr->z);
w121_2=R121*(z1-op_ptr->z);
w211_2=R211*(z1-op_ptr->z);
w111_2=R111*(z1-op_ptr->z);
G222=-Gamma*atan2(w222_1,w222_2)+Alpha*log(R222+y2-op_ptr->y)+Beta*log(R222+x2-op_ptr->x);
G122=-Gamma*atan2(w122_1,w122_2)+Alpha*log(R122+y2-op_ptr->y)+Beta*log(R122+x1-op_ptr->x);
G212=-Gamma*atan2(w212_1,w212_2)+Alpha*log(R212+y1-op_ptr->y)+Beta*log(R212+x2-op_ptr->x);
G112=-Gamma*atan2(w112_1,w112_2)+Alpha*log(R112+y1-op_ptr->y)+Beta*log(R112+x1-op_ptr->x);
G221=-Gamma*atan2(w221_1,w221_2)+Alpha*log(R221+y2-op_ptr->y)+Beta*log(R221+x2-op_ptr->x);
G121=-Gamma*atan2(w121_1,w121_2)+Alpha*log(R121+y2-op_ptr->y)+Beta*log(R121+x1-op_ptr->x);
G211=-Gamma*atan2(w211_1,w211_2)+Alpha*log(R211+y1-op_ptr->y)+Beta*log(R211+x2-op_ptr->x);
G111=-Gamma*atan2(w111_1,w111_2)+Alpha*log(R111+y1-op_ptr->y)+Beta*log(R111+x1-op_ptr->x);
return GCTL_T0/(4*GCTL_Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
}
Reference in New Issue View Git Blame Copy Permalink