From 8dffa168762847a44fdeb15e9322370a102013f3 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Wed, 5 Mar 2025 16:26:24 +0800 Subject: [PATCH] update --- lib/potential/gm_data.cpp | 104 +++++++++++++++++++++- lib/potential/gm_data.h | 39 +++++++- lib/potential/mkernel_tricone.cpp | 1 + lib/potential/mkernel_tricone_Ren2017.cpp | 3 +- 4 files changed, 142 insertions(+), 5 deletions(-) diff --git a/lib/potential/gm_data.cpp b/lib/potential/gm_data.cpp index 0afc35f..091c9b0 100644 --- a/lib/potential/gm_data.cpp +++ b/lib/potential/gm_data.cpp @@ -124,8 +124,8 @@ void gctl::read_IGRF_table(std::string file, array &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 &Mag_tensors, const return; } -void gctl::power_spectrum_shc(array &P, const array &C, const array &S, double R, double r, int n, int N) +void gctl::schmidt_legendre(double co_lati, int n_max, array &P, array &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 &P, const array &C, const array &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) { diff --git a/lib/potential/gm_data.h b/lib/potential/gm_data.h index 20d620f..c48ae07 100644 --- a/lib/potential/gm_data.h +++ b/lib/potential/gm_data.h @@ -198,6 +198,43 @@ namespace gctl */ void magnetic_tensors2deltaTs_sph(const array &Mag_tensors, const _1d_array &T0_inclina, const _1d_array &T0_declina, array &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 &P, array &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 &P, const array &C, const array &S, + void power_spectrum(array &P, const array &C, const array &S, double R, double r, int n, int N); } diff --git a/lib/potential/mkernel_tricone.cpp b/lib/potential/mkernel_tricone.cpp index 503e782..c6bf78b 100644 --- a/lib/potential/mkernel_tricone.cpp +++ b/lib/potential/mkernel_tricone.cpp @@ -25,6 +25,7 @@ * Also add information on how to contact you by electronic and paper mail. ******************************************************/ +#include #include "mkernel_tricone.h" void gctl::callink_magnetic_para(array &in_cone, array &out_para, const array &mag_B) diff --git a/lib/potential/mkernel_tricone_Ren2017.cpp b/lib/potential/mkernel_tricone_Ren2017.cpp index cc03f9d..bccc695 100644 --- a/lib/potential/mkernel_tricone_Ren2017.cpp +++ b/lib/potential/mkernel_tricone_Ren2017.cpp @@ -25,8 +25,9 @@ * Also add information on how to contact you by electronic and paper mail. ******************************************************/ +#include +#include #include "mkernel_tricone_Ren2017.h" -#include "cmath" double mag_tolerance = 1e-16;