zhangyiss/gctl_potential
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update

Browse Source
This commit is contained in:
张壹 2025-03-05 16:26:24 +08:00
parent b68d27f726
commit 8dffa16876
4 changed files with 142 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

104
lib/potential/gm_data.cpp
View File

@ -124,8 +124,8 @@ void gctl::read_IGRF_table(std::string file, array<IGRF_para> &IGRFs, int head_r
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;
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;
@ -241,7 +241,105 @@ void gctl::magnetic_tensors2deltaTs_sph(const array<tensor> &Mag_tensors, const
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)
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)
{

39
lib/potential/gm_data.h
View File

@ -198,6 +198,43 @@ namespace gctl
*/
void magnetic_tensors2deltaTs_sph(const array<tensor> &Mag_tensors, const _1d_array &T0_inclina, const _1d_array &T0_declina, array<point3dc> &deltaTs);
/**
* @brief This function evaluates all of the Schmidt-semi normalized associated Legendre
* functions up to degree nMax. The functions are initially scaled by
* 10^280 sin^m in order to minimize the effects of underflow at large m
* near the poles (see Holmes and Featherstone 2002, J. Geodesy, 76, 279-299).
* Note that this function performs the same operation as MAG_PcupLow.
* However this function also can be used for high degree (large nMax) models.
*
* Added into the GCTL by Yi Zhang in March 5th, 2025.
*
* Adopted from the FORTRAN code written by Mark Wieczorek September 25, 2005.
*
* Manoj Nair, Nov, 2009 Manoj.C.Nair@Noaa.Gov
*
* Change from the previous version
* The prevous version computes the derivatives as
* dP(n,m)(x)/dx, where x = sin(latitude) (or cos(colatitude) ).
* However, the WMM Geomagnetic routines requires dP(n,m)(x)/dlatitude.
* Hence the derivatives are multiplied by sin(latitude).
* Removed the options for CS phase and normalizations.
*
* Note: In geomagnetism, the derivatives of ALF are usually found with
* respect to the colatitudes. Here the derivatives are found with respect
* to the latitude. The difference is a sign reversal for the derivative of
* the Associated Legendre Functions.
*
* The derivatives can't be computed for latitude = |90| degrees.
*
* @param co_lati Co-latitude in degrees.
* @param n_max Maximum spherical harmonic degree to compute.
* @param P A vector of all associated Legendgre polynomials evaluated at co_lati up to nMax.
* The lenght must by greater or equal to (n_max+1)*(n_max+2)/2. The coefficients are stored
* in a degree major order.
* @param dP Derivative of P with respect to latitude
*/
void schmidt_legendre(double co_lati, int n_max, array<double> &P, array<double> &dP);
/**
* @brief Calculate the power spectrum of given spherical harmonic coefficients.
*
@ -209,7 +246,7 @@ namespace gctl
* @param n Strat order
* @param N End order
*/
void power_spectrum_shc(array<double> &P, const array<double> &C, const array<double> &S,
void power_spectrum(array<double> &P, const array<double> &C, const array<double> &S,
double R, double r, int n, int N);
}

1
lib/potential/mkernel_tricone.cpp
View File

@ -25,6 +25,7 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include <map>
#include "mkernel_tricone.h"
void gctl::callink_magnetic_para(array<magcone> &in_cone, array<magcone_para> &out_para, const array<point3dc> &mag_B)

3
lib/potential/mkernel_tricone_Ren2017.cpp
View File

@ -25,8 +25,9 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include <cmath>
#include <map>
#include "mkernel_tricone_Ren2017.h"
#include "cmath"
double mag_tolerance = 1e-16;