265 lines
12 KiB
C++
265 lines
12 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 "gm_data.h"
|
|
|
|
gctl::tensor gctl::transform_matrix(const point3ds &op)
|
|
{
|
|
tensor R;
|
|
|
|
R[0][0] = sin((0.5 - op.lat/180.0)*GCTL_Pi)*cos((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[0][1] = sin((0.5 - op.lat/180.0)*GCTL_Pi)*sin((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[0][2] = cos((0.5 - op.lat/180.0)*GCTL_Pi);
|
|
|
|
R[1][0] = cos((0.5 - op.lat/180.0)*GCTL_Pi)*cos((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[1][1] = cos((0.5 - op.lat/180.0)*GCTL_Pi)*sin((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[1][2] = -1.0*sin((0.5 - op.lat/180.0)*GCTL_Pi);
|
|
|
|
R[2][0] = -1.0*sin((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[2][1] = cos((2.0 + op.lon/180.0)*GCTL_Pi);
|
|
R[2][2] = 0.0;
|
|
|
|
return R;
|
|
}
|
|
|
|
void gctl::geomag_local2Cartesian(array<point3dc> &abs_b_, const array<point3dc> &geo_b_, const array<point3ds> &loc_s_)
|
|
{
|
|
if (geo_b_.size() != loc_s_.size())
|
|
{
|
|
throw std::runtime_error("[gctl::geomag_local2Cartesian] Incompatible input arraies' size.");
|
|
}
|
|
|
|
abs_b_.resize(geo_b_.size());
|
|
|
|
point3dc mag_v;
|
|
for (size_t i = 0; i < geo_b_.size(); i++)
|
|
{
|
|
// magnetization vector at the local coordinates
|
|
// note here the postive direction of the inclination angle is downward
|
|
// note here the postive direction of the declination angle is clockwise
|
|
mag_v.x = geo_b_[i].z;
|
|
mag_v.y = geo_b_[i].y;
|
|
mag_v.z = -1.0*geo_b_[i].x;
|
|
|
|
// magnetic susceptibility is taken as one here
|
|
// rotate the local coordinate system to the regular status
|
|
abs_b_[i] = mag_v.rotate((90.0 - loc_s_[i].lat)*GCTL_Pi/180.0, 0.0, (90.0 + loc_s_[i].lon)*GCTL_Pi/180.0);
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::read_IGRF_table(std::string file, array<IGRF_para> &IGRFs, int head_record, std::string ext_f)
|
|
{
|
|
std::ifstream infile;
|
|
gctl::open_infile(infile, file, ext_f);
|
|
|
|
std::vector<IGRF_para> vec_para;
|
|
IGRF_para tmp_para;
|
|
|
|
std::string err_str, tmp_str;
|
|
std::stringstream tmp_ss;
|
|
std::string date_str, alti_str, dd_str, dm_str, id_str, im_str;
|
|
|
|
// 头信息就跳过
|
|
for (size_t i = 0; i < head_record; i++)
|
|
{
|
|
getline(infile, tmp_str);
|
|
}
|
|
|
|
while (getline(infile, tmp_str))
|
|
{
|
|
//Date Coord-System Altitude Latitude Longitude D_deg D_min I_deg I_min H_nT X_nT Y_nT Z_nT F_nT dD_min dI_min dH_nT dX_nT dY_nT dZ_nT dF_nT
|
|
//2017,6,15 D M4000.00 30.0167 105.017 -2d 16m 46d 52m 34358.9 34332.2 -1354.1 36667.2 50249.6 -3.8 5.6 -13.2 -14.8 -37.9 105.0 67.6
|
|
tmp_ss.clear();
|
|
tmp_ss.str(tmp_str);
|
|
tmp_ss >> date_str >> tmp_para.CSystem >> alti_str
|
|
>> tmp_para.Latitude >> tmp_para.Longitude
|
|
>> dd_str >> dm_str >> id_str >> im_str
|
|
>> tmp_para.H_nT >> tmp_para.X_nT >> tmp_para.Y_nT
|
|
>> tmp_para.Z_nT >> tmp_para.F_nT >> tmp_para.dD_min
|
|
>> tmp_para.dI_min >> tmp_para.dH_nT >> tmp_para.dX_nT
|
|
>> tmp_para.dY_nT >> tmp_para.dZ_nT >> tmp_para.dF_nT;
|
|
|
|
if (1 != sscanf(dd_str.c_str(), "%dd", &tmp_para.D_deg) ||
|
|
1 != sscanf(dm_str.c_str(), "%dm", &tmp_para.D_min) ||
|
|
1 != sscanf(id_str.c_str(), "%dd", &tmp_para.I_deg) ||
|
|
1 != sscanf(im_str.c_str(), "%dm", &tmp_para.I_min) ||
|
|
2 != sscanf(alti_str.c_str(), "%c%lf", &tmp_para.AltiCode, &tmp_para.Altitude))
|
|
{
|
|
throw gctl::invalid_argument("[gctl::read_IGRF_table] Wrong format: " + tmp_str);
|
|
}
|
|
|
|
if (3 != sscanf(date_str.c_str(), "%d,%d,%d", &tmp_para.Year, &tmp_para.Month, &tmp_para.Day))
|
|
{
|
|
double deci_year;
|
|
if (1 != sscanf(date_str.c_str(), "%lf", &deci_year))
|
|
{
|
|
throw gctl::invalid_argument("[gctl::read_IGRF_table] Wrong format: " + tmp_str);
|
|
}
|
|
|
|
decimal_year2date(deci_year, tmp_para.Year, tmp_para.Month, tmp_para.Day);
|
|
}
|
|
|
|
tmp_para.I = fabs(tmp_para.I_deg) + tmp_para.I_min/60.0;
|
|
tmp_para.D = fabs(tmp_para.D_deg) + tmp_para.D_min/60.0;
|
|
if (tmp_para.I_deg < 0) tmp_para.I *= -1;
|
|
if (tmp_para.D_deg < 0) tmp_para.D *= -1;
|
|
|
|
vec_para.push_back(tmp_para);
|
|
}
|
|
infile.close();
|
|
|
|
IGRFs.input(vec_para);
|
|
destroy_vector(vec_para);
|
|
return;
|
|
}
|
|
|
|
void gctl::magnetic_components2deltaT(const _1d_array &Hax, const _1d_array &Hay, const _1d_array &Za,
|
|
_1d_array &deltaT, double T0_inclina, double T0_declina)
|
|
{
|
|
if (Hax.size() != Hay.size() || Hay.size() != Za.size())
|
|
{
|
|
throw std::runtime_error("[gctl::magnetic_components2deltaT] Unmatched components' size.");
|
|
}
|
|
|
|
// note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
|
|
double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina/180.0);
|
|
|
|
deltaT.resize(Hax.size());
|
|
for (size_t i = 0; i < Hax.size(); i++)
|
|
{
|
|
deltaT[i] = hax_f*Hax[i] + hay_f*Hay[i] + za_f*Za[i];
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::magnetic_components2deltaT(const array<point3dc> &Mag_components, _1d_array &deltaT, double T0_inclina, double T0_declina)
|
|
{
|
|
// note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
|
|
double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina/180.0);
|
|
|
|
deltaT.resize(Mag_components.size());
|
|
for (size_t i = 0; i < Mag_components.size(); i++)
|
|
{
|
|
deltaT[i] = hax_f*Mag_components[i].x + hay_f*Mag_components[i].y + za_f*Mag_components[i].z;
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::magnetic_tensors2deltaTs(const array<tensor> &Mag_tensors, array<point3dc> &deltaTs, double T0_inclina, double T0_declina)
|
|
{
|
|
// note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
|
|
double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina/180.0);
|
|
|
|
deltaTs.resize(Mag_tensors.size());
|
|
for (size_t i = 0; i < Mag_tensors.size(); i++)
|
|
{
|
|
deltaTs[i].x = hax_f*Mag_tensors[i][0][0] + hay_f*Mag_tensors[i][0][1] + za_f*Mag_tensors[i][0][2];
|
|
deltaTs[i].y = hax_f*Mag_tensors[i][0][1] + hay_f*Mag_tensors[i][1][1] + za_f*Mag_tensors[i][1][2];
|
|
deltaTs[i].z = hax_f*Mag_tensors[i][0][2] + hay_f*Mag_tensors[i][1][2] + za_f*Mag_tensors[i][2][2];
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::magnetic_components2deltaT_sph(const array<point3dc> &Mag_components, const _1d_array &T0_inclina, const _1d_array &T0_declina, _1d_array &deltaT)
|
|
{
|
|
// note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
|
|
deltaT.resize(Mag_components.size());
|
|
for (size_t i = 0; i < Mag_components.size(); i++)
|
|
{
|
|
double hax_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*sin(GCTL_Pi*T0_declina[i]/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*cos(GCTL_Pi*T0_declina[i]/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina[i]/180.0);
|
|
|
|
// note here Mag_components[i].x is the reversed radial component
|
|
// Mag_components[i].y is the latitudinal component (south pointing). to use it in a local cratesian coordinate, we need to reverse it
|
|
// Mag_components[i].z is the longtidinal component (east pointing)
|
|
deltaT[i] = hax_f*Mag_components[i].z - hay_f*Mag_components[i].y + za_f*Mag_components[i].x;
|
|
}
|
|
return;
|
|
}
|
|
|
|
double gctl::magnetic_components2deltaT_sph(const point3dc &Mag_components, double T0_inclina, double T0_declina)
|
|
{
|
|
double hax_f = cos(GCTL_Pi*T0_inclina/180.0)*sin(GCTL_Pi*T0_declina/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina/180.0)*cos(GCTL_Pi*T0_declina/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina/180.0);
|
|
|
|
// note here Mag_components[i].x is the reversed radial component
|
|
// Mag_components[i].y is the latitudinal component (south pointing). to use it in a local cratesian coordinate, we need to reverse it
|
|
// Mag_components[i].z is the longtidinal component (east pointing)
|
|
return hax_f*Mag_components.z - hay_f*Mag_components.y + za_f*Mag_components.x;
|
|
}
|
|
|
|
void gctl::magnetic_tensors2deltaTs_sph(const array<tensor> &Mag_tensors, const _1d_array &T0_inclina, const _1d_array &T0_declina, array<point3dc> &deltaTs)
|
|
{
|
|
// note that we use regular right-hand cartesian coordinates here, not the reversed z-axis one.
|
|
deltaTs.resize(Mag_tensors.size());
|
|
for (size_t i = 0; i < Mag_tensors.size(); i++)
|
|
{
|
|
double hax_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*sin(GCTL_Pi*T0_declina[i]/180.0);
|
|
double hay_f = cos(GCTL_Pi*T0_inclina[i]/180.0)*cos(GCTL_Pi*T0_declina[i]/180.0);
|
|
double za_f = sin(GCTL_Pi*T0_inclina[i]/180.0);
|
|
|
|
// note here Mag_tensors[i][:][0] is the reversed radial component
|
|
// Mag_components[i][:][1] is the latitudinal component (south pointing). to use it in a local cratesian coordinate, we need to reverse it
|
|
// Mag_components[i][:][2] is the longtidinal component (east pointing)
|
|
deltaTs[i].x = za_f*Mag_tensors[i][0][0] - hay_f*Mag_tensors[i][0][1] + hax_f*Mag_tensors[i][0][2];
|
|
deltaTs[i].y = za_f*Mag_tensors[i][0][1] - hay_f*Mag_tensors[i][1][1] + hax_f*Mag_tensors[i][1][2];
|
|
deltaTs[i].z = za_f*Mag_tensors[i][0][2] - hay_f*Mag_tensors[i][1][2] + hax_f*Mag_tensors[i][2][2];
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::power_spectrum_shc(array<double> &P, const array<double> &C, const array<double> &S, double R, double r, int n, int N)
|
|
{
|
|
if (n < 0 || n >= N || C.size() < (N + 1)*(N + 2)/2 || S.size() < (N + 1)*(N + 2)/2)
|
|
{
|
|
throw std::runtime_error("[gctl::power_spectrum_shc] Invalid input parameters.");
|
|
}
|
|
|
|
P.resize(N - n + 1, 0);
|
|
|
|
int id;
|
|
for (int i = n; i < N + 1; i++)
|
|
{
|
|
id = i*(i + 1)/2;
|
|
for (int j = 0; j < i + 1; j++)
|
|
{
|
|
P[i - n] += (C[id + j]*C[id + j] + S[id + j]*S[id + j]);
|
|
}
|
|
|
|
P[i - n] *= (n + 1)*pow(R/r, 2*n + 4);
|
|
}
|
|
return;
|
|
} |