469 lines
17 KiB
C++
469 lines
17 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 = abs(tmp_para.I_deg) + tmp_para.I_min/60.0;
|
|
tmp_para.D = abs(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::read_Swarm_shc(std::string file, array<shc_data> &SHCs, int &spline_odr, int &sample_num)
|
|
{
|
|
dsv_io tio;
|
|
tio.head_number(2);
|
|
tio.load_text(file);
|
|
|
|
// parse the head records
|
|
int n, N, tnum;
|
|
std::vector<double> time_snaps;
|
|
std::vector<std::string> rcd = tio.head_records();
|
|
parse_string_to_value(rcd[0], ' ', true, n, N, tnum, spline_odr, sample_num);
|
|
parse_string_to_vector(rcd[1], ' ', time_snaps);
|
|
|
|
array<int> nid, mid;
|
|
array<double> vals;
|
|
int rnum = tio.row_number();
|
|
tio.get_column(nid, 1);
|
|
tio.get_column(mid, 2);
|
|
|
|
// store coefficients
|
|
int id, shc_num;
|
|
SHCs.resize(tnum);
|
|
for (size_t i = 0; i < tnum; i++)
|
|
{
|
|
SHCs[i].ns = n;
|
|
SHCs[i].ms = 0;
|
|
SHCs[i].ne = N;
|
|
SHCs[i].me = N;
|
|
SHCs[i].time = time_snaps[i];
|
|
|
|
shc_num = (N + 1)*(N + 2)/2;
|
|
SHCs[i].Snm.resize(shc_num, 0.0);
|
|
SHCs[i].Cnm.resize(shc_num, 0.0);
|
|
|
|
tio.get_column(vals, i + 3);
|
|
for (size_t j = 0; j < rnum; j++)
|
|
{
|
|
if (mid[j] > 0)
|
|
{
|
|
id = nid[j]*(nid[j] + 1)/2 + mid[j];
|
|
SHCs[i].Snm[id] = vals[j];
|
|
SHCs[i].Cnm[id] = vals[j + 1];
|
|
}
|
|
else if (mid[j] == 0)
|
|
{
|
|
id = nid[j]*(nid[j] + 1)/2;
|
|
SHCs[i].Snm[id] = vals[j];
|
|
SHCs[i].Cnm[id] = 0.0;
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::save_Swarm_shc(std::string file, const array<shc_data> &SHCs, int spline_odr, int sample_num, std::string info)
|
|
{
|
|
time_t now = time(0);
|
|
char* dt = ctime(&now);
|
|
|
|
std::ofstream ofile;
|
|
open_outfile(ofile, file + ".shc", ".txt");
|
|
ofile << "# " << info << std::endl;
|
|
ofile << "# Generated on " << dt;
|
|
|
|
int s = SHCs[0].ns, S = SHCs[0].ne;
|
|
int c = SHCs.size();
|
|
ofile << s << " " << S << " " << c << " " << spline_odr << " " << sample_num << std::endl;
|
|
|
|
ofile << SHCs[0].time;
|
|
for (size_t i = 1; i < SHCs.size(); i++)
|
|
{
|
|
ofile << " " << SHCs[i].time;
|
|
}
|
|
|
|
for (size_t n = s; n <= S; n++)
|
|
{
|
|
for (size_t m = 0; m <= n; m++)
|
|
{
|
|
if (m == 0)
|
|
{
|
|
ofile << "\n" << n << " " << m;
|
|
for (size_t i = 0; i < c; i++)
|
|
{
|
|
ofile << " " << SHCs[i].coeff_s(n, m);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
ofile << "\n" << n << " " << m;
|
|
for (size_t i = 0; i < c; i++)
|
|
{
|
|
ofile << " " << SHCs[i].coeff_s(n, m);
|
|
}
|
|
|
|
ofile << "\n" << n << " -" << m;
|
|
for (size_t i = 0; i < c; i++)
|
|
{
|
|
ofile << " " << SHCs[i].coeff_c(n, m);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
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::schmidt_legendre(double co_lati, int n_max, array<double> &P, array<double> &dP)
|
|
{
|
|
if (n_max <= 0) throw std::runtime_error("[gctl::schmidt_legendre] Invalid paramters.");
|
|
|
|
double cos_lat = cosd(co_lati);
|
|
if (abs(cos_lat) == 1.0) throw std::runtime_error("[gctl::schmidt_legendre] Derivative cannot be calculated at poles.");
|
|
|
|
double pm2, pm1, pmm, plm, rescalem, z, scalef;
|
|
int k, kstart, m, n, NumTerms;
|
|
|
|
NumTerms = ((n_max + 1) * (n_max + 2) / 2);
|
|
P.resize(NumTerms, 0.0);
|
|
dP.resize(NumTerms, 0.0);
|
|
_1d_array f1(NumTerms + 1), f2(NumTerms + 1), PreSqr(NumTerms + 1);
|
|
|
|
scalef = 1.0e-280;
|
|
|
|
for(n = 0; n <= 2 * n_max + 1; ++n)
|
|
{
|
|
PreSqr[n] = sqrt((double) (n));
|
|
}
|
|
|
|
k = 2;
|
|
for(n = 2; n <= n_max; n++)
|
|
{
|
|
k = k + 1;
|
|
f1[k] = (double) (2 * n - 1) / (double) (n);
|
|
f2[k] = (double) (n - 1) / (double) (n);
|
|
for(m = 1; m <= n - 2; m++)
|
|
{
|
|
k = k + 1;
|
|
f1[k] = (double) (2 * n - 1) / PreSqr[n + m] / PreSqr[n - m];
|
|
f2[k] = PreSqr[n - m - 1] * PreSqr[n + m - 1] / PreSqr[n + m] / PreSqr[n - m];
|
|
}
|
|
k = k + 2;
|
|
}
|
|
|
|
/*z = sin (geocentric latitude) */
|
|
z = sqrt((1.0 - cos_lat)*(1.0 + cos_lat));
|
|
pm2 = 1.0;
|
|
P[0] = 1.0;
|
|
dP[0] = 0.0;
|
|
|
|
pm1 = cos_lat;
|
|
P[1] = pm1;
|
|
dP[1] = z;
|
|
k = 1;
|
|
|
|
for(n = 2; n <= n_max; n++)
|
|
{
|
|
k = k + n;
|
|
plm = f1[k] * cos_lat * pm1 - f2[k] * pm2;
|
|
P[k] = plm;
|
|
dP[k] = (double) (n) * (pm1 - cos_lat * plm) / z;
|
|
pm2 = pm1;
|
|
pm1 = plm;
|
|
}
|
|
|
|
pmm = PreSqr[2] * scalef;
|
|
rescalem = 1.0 / scalef;
|
|
kstart = 0;
|
|
|
|
for(m = 1; m <= n_max - 1; ++m)
|
|
{
|
|
rescalem = rescalem*z;
|
|
|
|
/* Calculate Pcup(m,m)*/
|
|
kstart = kstart + m + 1;
|
|
pmm = pmm * PreSqr[2 * m + 1] / PreSqr[2 * m];
|
|
P[kstart] = pmm * rescalem / PreSqr[2 * m + 1];
|
|
dP[kstart] = -((double) (m) * cos_lat * P[kstart] / z);
|
|
pm2 = pmm / PreSqr[2 * m + 1];
|
|
/* Calculate Pcup(m+1,m)*/
|
|
k = kstart + m + 1;
|
|
pm1 = cos_lat * PreSqr[2 * m + 1] * pm2;
|
|
P[k] = pm1*rescalem;
|
|
dP[k] = ((pm2 * rescalem) * PreSqr[2 * m + 1] - cos_lat * (double) (m + 1) * P[k]) / z;
|
|
/* Calculate Pcup(n,m)*/
|
|
for(n = m + 2; n <= n_max; ++n)
|
|
{
|
|
k = k + n;
|
|
plm = cos_lat * f1[k] * pm1 - f2[k] * pm2;
|
|
P[k] = plm*rescalem;
|
|
dP[k] = (PreSqr[n + m] * PreSqr[n - m] * (pm1 * rescalem) - (double) (n) * cos_lat * P[k]) / z;
|
|
pm2 = pm1;
|
|
pm1 = plm;
|
|
}
|
|
}
|
|
|
|
/* Calculate Pcup(nMax,nMax)*/
|
|
rescalem = rescalem*z;
|
|
kstart = kstart + m + 1;
|
|
pmm = pmm / PreSqr[2 * n_max];
|
|
P[kstart] = pmm * rescalem;
|
|
dP[kstart] = -(double) (n_max) * cos_lat * P[kstart] / z;
|
|
return;
|
|
}
|
|
|
|
void gctl::power_spectrum(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] 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] *= (i + 1)*pow(R/r, 2*i + 4);
|
|
}
|
|
return;
|
|
} |