/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * Geophysical Computational Tools & Library (GCTL) * * Copyright (c) 2023 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 . * * 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 "legendre.h" double gctl::legendre_polynomials(size_t order, double x, bool derivative) { //if (x > 1.0 || x < -1.0) //{ // throw gctl::runtime_error("The calculating position must be in [-1, 1]. From gctl::legendre_polynomials(...)"); //} if (derivative == false) { if (order == 0) return 1.0; if (order == 1) return x; double n = 0, p = 0, p_1 = x, p_2 = 1.0; for (size_t i = 2; i <= order; i++) { n = (double) i - 1.0; p = (2.0*n + 1.0)*x*p_1/(n + 1.0) - n*p_2/(n + 1.0); p_2 = p_1; p_1 = p; } return p; } else { if (order == 0) return 0.0; if (order == 1) return 1.0; double n = 0, p = 0, p_1 = x, p_2 = 1.0; double pd = 0, pd_1 = 1.0, pd_2 = 0.0; for (size_t i = 2; i <= order; i++) { n = (double) i - 1.0; p = (2.0*n + 1.0)*x*p_1/(n + 1.0) - n*p_2/(n + 1.0); pd = (2.0*n + 1.0)*(p_1 + x*pd_1)/(n + 1.0) - n*pd_2/(n + 1.0); p_2 = p_1; p_1 = p; pd_2 = pd_1; pd_1 = pd; } return pd; } } void gctl::get_a_nm_array(int max_order, array> &cs) { int i, j; if (cs.size() != max_order) { cs.resize(max_order); for (i = 0; i < max_order; i++) cs[i].resize(i+1); } //向下列推计算 #pragma omp parallel for private(i,j) schedule(guided) for (j = 0; j < max_order; j++) { cs[j][j] = 0; //对角线上的值直接给0 反正用不到 for (i = j+1; i < max_order; i++) { cs[i][j] = sqrt(((2.0*i-1)*(2.0*i+1))/((i-j)*(i+j))); } } return; } void gctl::get_b_nm_array(int max_order, array> &cs) { int i,j; if (cs.size() != max_order) { cs.resize(max_order); for (i = 0; i < max_order; i++) cs[i].resize(i+1); } //向下列推计算 #pragma omp parallel for private(i,j) schedule(guided) for (j = 0; j < max_order; j++) { cs[j][j] = 0; //对角线上的值直接给0 反正用不到 for (i = j+1; i < max_order; i++) { cs[i][j] = sqrt(((2.0*i+1)*(i+j-1)*(i-j-1))/((i-j)*(i+j)*(2.0*i-3))); } } return; } void gctl::nalf_sfcm(array> &nalf, const array> &a_nm, const array> &b_nm, int max_order, double theta, legendre_norm_e norm, bool derivative) { if (a_nm.size() != max_order || b_nm.size() != max_order) { std::string err_str = "Incompatible coefficients' size."; throw runtime_error(err_str); } double norSum; if (norm == One) { norSum = 1.0; } else if (norm == Pi4) { norSum = 4.0*GCTL_Pi; } if (nalf.size() != max_order) { nalf.resize(max_order); for (int i = 0; i < max_order; i++) nalf[i].resize(i+1); } if (derivative) { double tmp; array tmp_nalf(max_order, 0.0), tmp2_nalf(max_order, 0.0); nalf[0][0] = 0.0; tmp_nalf[0] = sqrt(norSum)/sqrt(4.0*GCTL_Pi); nalf[1][1] = sqrt(3.0)*cos(theta*GCTL_Pi/180.0); tmp_nalf[1] = sqrt(3.0)*sin(theta*GCTL_Pi/180.0); //计算对角线上的值递归计算不能并行 for (int i = 2; i < max_order; i++) { nalf[i][i] = cos(theta*GCTL_Pi/180.0)*sqrt(0.5*(2.0*i+1)/i)*tmp_nalf[i-1] + sin(theta*GCTL_Pi/180.0)*sqrt(0.5*(2.0*i+1)/i)*nalf[i-1][i-1]; tmp_nalf[i] = sin(theta*GCTL_Pi/180.0)*sqrt(0.5*(2.0*i+1)/i)*tmp_nalf[i-1]; } //计算次对角线(m+1,m)上的值递归计算不能并行 for (int i = 0; i < max_order-1; i++) { nalf[i+1][i] = -1.0*sin(theta*GCTL_Pi/180.0)*sqrt(2.0*i+3)*tmp_nalf[i] + cos(theta*GCTL_Pi/180.0)*sqrt(2.0*i+3)*nalf[i][i]; tmp2_nalf[i] = tmp_nalf[i]; tmp_nalf[i] = cos(theta*GCTL_Pi/180.0)*sqrt(2.0*i+3)*tmp_nalf[i]; } //这里可以使用并行加速计算外层循环内层计算因为是递归计算因此不能并行 int i,j; #pragma omp parallel for private(i,j, tmp) schedule(guided) for (j = 0; j < max_order-1; j++) { for (i = j+2; i < max_order; i++) { nalf[i][j] = -1.0*a_nm[i][j]*sin(theta*GCTL_Pi/180.0)*tmp_nalf[j] + a_nm[i][j]*cos(theta*GCTL_Pi/180.0)*nalf[i-1][j] - b_nm[i][j]*nalf[i-2][j]; tmp = tmp_nalf[j]; tmp_nalf[j] = a_nm[i][j]*cos(theta*GCTL_Pi/180.0)*tmp_nalf[j] - b_nm[i][j]*tmp2_nalf[j]; tmp2_nalf[j] = tmp; } } } else { //赋初值给前两个对角线上的值 //norSum为1时第一个值为sqrt(1)/sqrt(4.0*pi) = 1/sqrt(4.0*pi), 对应的归一化值为1 //norSum为4.0*pi时第一个值为4.0*pi/sqrt(4.0*pi) = 1, 对应的归一化值为4.0*pi nalf[0][0] = sqrt(norSum)/sqrt(4.0*GCTL_Pi); nalf[1][1] = sqrt(3.0)*sin(theta*GCTL_Pi/180.0); //计算对角线上的值递归计算不能并行 for (int i = 2; i < max_order; i++) { nalf[i][i] = sin(theta*GCTL_Pi/180.0)*sqrt(0.5*(2.0*i+1)/i)*nalf[i-1][i-1]; } //计算次对角线(m+1,m)上的值递归计算不能并行 for (int i = 0; i < max_order-1; i++) { nalf[i+1][i] = cos(theta*GCTL_Pi/180.0)*sqrt(2.0*i+3)*nalf[i][i]; } //这里可以使用并行加速计算外层循环内层计算因为是递归计算因此不能并行 int i,j; #pragma omp parallel for private(i,j) schedule(guided) for (j = 0; j < max_order-1; j++) { for (i = j+2; i < max_order; i++) { nalf[i][j] = a_nm[i][j]*cos(theta*GCTL_Pi/180.0)*nalf[i-1][j] - b_nm[i][j]*nalf[i-2][j]; } } } return; }