112 lines
4.1 KiB
C++
112 lines
4.1 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 "mkernel_dipole.h"
|
|
|
|
void gctl::mag_axial_dipole_parameters(double &M, double &I, double r, double lat_deg, double sus, double ref_r, double abs_g0)
|
|
{
|
|
M = sus*fabs(abs_g0)*power3(ref_r/r)*sqrt(1.0 + 3.0*power2(sind(lat_deg)));
|
|
I = atan(2.0*tan(GCTL_Pi*lat_deg/180.0));
|
|
return;
|
|
}
|
|
|
|
gctl::point3dc gctl::magkernel_single(const mag_dipole &a_dipole, const point3dc &a_op)
|
|
{
|
|
double r = a_op.module();
|
|
// T -> nT 1e+9 * 1e-7 -> 1e+2
|
|
return 1e+2*a_dipole.M*(3.0*dot(a_dipole.n, a_op)/power5(r)*a_op - 1.0/power3(r)*a_dipole.n);
|
|
}
|
|
|
|
gctl::point3dc gctl::magkernel_single(const mag_dipole &a_dipole, const point3ds &a_op, tensor *R_ptr)
|
|
{
|
|
tensor R;
|
|
if (R_ptr != nullptr) R = *R_ptr;
|
|
else R = transform_matrix(a_op);
|
|
|
|
double r = a_op.rad;
|
|
point3dc pc = a_op.s2c();
|
|
|
|
// T -> nT 1e+9 * 1e-7 -> 1e+2
|
|
point3dc b = 1e+2*a_dipole.M*(3.0*dot(a_dipole.n, pc)/power5(r)*pc - 1.0/power3(r)*a_dipole.n);
|
|
b = R*b;
|
|
b.z = 0.0; // the latitudinal componement is zero
|
|
return b;
|
|
}
|
|
|
|
void gctl::magobser(array<point3dc> &out_obs, const array<mag_dipole> &dipoles,
|
|
const array<point3dc> &obsp, verbose_type_e verbose)
|
|
{
|
|
int i, j;
|
|
int o_size = obsp.size();
|
|
int e_size = dipoles.size();
|
|
|
|
out_obs.resize(o_size, point3dc(0.0, 0.0, 0.0));
|
|
|
|
gctl::progress_bar bar(e_size, "magobser_componments");
|
|
for (j = 0; j < e_size; j++)
|
|
{
|
|
if (verbose == gctl::FullMsg) bar.progressed(j);
|
|
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
|
|
|
#pragma omp parallel for private (i) schedule(guided)
|
|
for (i = 0; i < o_size; i++)
|
|
{
|
|
out_obs[i] = out_obs[i] + magkernel_single(dipoles[j], obsp[i]);;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::magobser(array<point3dc> &out_obs, const array<mag_dipole> &dipoles,
|
|
const array<point3ds> &obsp, verbose_type_e verbose)
|
|
{
|
|
int i, j;
|
|
int o_size = obsp.size();
|
|
int e_size = dipoles.size();
|
|
|
|
out_obs.resize(o_size, point3dc(0.0, 0.0, 0.0));
|
|
array<tensor> Rs(o_size);
|
|
#pragma omp parallel for private (i) schedule(guided)
|
|
for (i = 0; i < o_size; i++)
|
|
{
|
|
Rs[i] = transform_matrix(obsp[i]);
|
|
}
|
|
|
|
gctl::progress_bar bar(e_size, "magobser_componments");
|
|
for (j = 0; j < e_size; j++)
|
|
{
|
|
if (verbose == gctl::FullMsg) bar.progressed(j);
|
|
else if (verbose == gctl::ShortMsg) bar.progressed_simple(j);
|
|
|
|
#pragma omp parallel for private (i) schedule(guided)
|
|
for (i = 0; i < o_size; i++)
|
|
{
|
|
out_obs[i] = out_obs[i] + magkernel_single(dipoles[j], obsp[i], Rs.get(i));
|
|
}
|
|
}
|
|
return;
|
|
} |