В библиотеке для научных и технических расчетов SciPy для языка Python есть возможность работы со многими специальными функциями, в том числе и со сферическими, которые реализованы в scipy.special.sph_harm(k, n, lon, colat) (обозначения изменены по сравнению с официальной документацией на более привычные). Здесь k - порядок (order), n - степень (degree), lon - долгота $\lambda$, colat - полярное расстояние $\theta$, то есть дополнение широты $\phi$ до $90^\circ$.

Сферическая функция степени $n$ и порядка $k$ определена в scipy.special.sph_harm довольно непривычным для геодезистов образом $$ Y_{n}^{k} (\theta, \lambda) = \sqrt{\frac{2n + 1}{4\pi}\frac{(n - k)!}{(n + k)!}} e^{ik\lambda} (-1)^k P_n^k (\cos{\theta}), \quad 0 \leq n < \infty, \quad -n \leq k \leq n. $$ Это нормированная комплексная сферическая функция с фазой Кондона-Шортли. Здесь $P_n^k (\cos{\theta})$ - присоединённая функция Лежандра, $i$ - мнимая единица. Теперь разберемся с каждым элементом.

Комплексная сферическая функция

Всякую вещественную квадратично интегрируемую функцию, заданную на сфере, можно представить в виде ряда комплексных сферических функций $$ f \left(\theta, \lambda\right) = \sum\limits_{n=0}^{\infty}\sum\limits_{k=-n}^{n} f_n^k Y_n^k, \quad Y_{n}^{k} \left(\theta, \lambda\right) = e^{ik\lambda} P_n^k (\cos{\theta}), $$ где $f_n^k$ - комплексные гармонические коэффициенты. Этот ряд приведем к более привычному вещественному виду $$ f (\theta, \lambda) = \sum\limits_{n=0}^{\infty}\sum\limits_{k=0}^{n} \left(a_{nk}\cos{k\lambda} + b_{nk}\sin{k\lambda}\right) P_n^k(\cos\theta), $$ где $a_{nk}, b_{nk}$ - вещественные гармонические коэффициенты.

Разобъём комплексный ряд на три части для $k < 0$, $k = 0$, $k > 0$, получим $$ \begin{align*} \sum\limits_{n=0}^{\infty}\sum\limits_{k=-n}^{n} f_n^k Y_n^k &= \sum\limits_{n=0}^{\infty}\left[ f_n^0 Y_n^0 + \sum\limits_{k=-n}^{-1} f_n^k Y_n^k + \sum\limits_{k=1}^{n} f_n^k Y_n^k \right] = \\\\ &= \sum\limits_{n=0}^{\infty}\left[ f_n^0 P_n + \sum\limits_{k=1}^{n} f_n^{-k} Y_n^{-k} + \sum\limits_{k=1}^{n} f_n^k Y_n^k \right]. \end{align*} $$ Затем воспользуемся формулой Эйлера $$ e^{ix} = \cos{x} - i\sin{x} $$ и получим $$ \begin{align*} &\sum\limits_{n=0}^{\infty}\left[ f_n^0 P_n + \sum\limits_{k=1}^{n} f_n^{-k} (\cos{k\lambda} + i\sin{k\lambda}) P_n^{-k} + \sum\limits_{k=1}^{n} f_n^k (\cos{k\lambda} - i\sin{k\lambda}) P_n^k \right] = \\\\ &\sum\limits_{n=0}^{\infty}\left[ f_n^0 P_n + \sum\limits_{k=1}^{n} \left[ \left( f_n^k + f_n^{-k} \right)\cos{k\lambda} + i\left( f_n^k - f_n^{-k} \right)\sin{k\lambda}\right] P_n^k \right]. \end{align*} $$

Введём следующие обозначения $$ a_{nk} = f_n^k + f_n^{-k},\qquad a_{n0} = 2f_n^0, \qquad b_{nk} = i\left( f_n^k - f_n^{-k} \right) $$ и окончательно получим, как и хотели, вещественный ряд Лапласа $$ f (\theta, \lambda) = \sum\limits_{n=0}^{\infty}\sum\limits_{k=0}^{n} (a_{nk}\cos{k\lambda} + b_{nk}\sin{k\lambda}) P_n^k(\cos\theta). $$

Его также можно рассматривать как вещественную часть ряда комплексных сферических функций с коэффициентами $f_n^k = a_{nk}- ib_{nk}$: $$ f (\theta, \lambda) = \sum\limits_{n=0}^{\infty}\sum\limits_{k=-n}^{n} f_n^k Y_n^k = \sum\limits_{n=0}^{\infty}\sum\limits_{k=-n}^{n} (a_{nk}- ib_{nk}) e^{ik\lambda} P_n^k(\cos\theta), $$ что легко показать снова с помощью формулы Эйлера.

Таким образом, если коэффициенты комплексных гармоник представлены через коэффициенты вещественных гармоник в виде $f_n^k = a_{nk} - ib_{nk}$, то вещественная часть произведения $f_n^k Y_n^k$ будет равна $$ f_n^k Y_n^k = \left(a_{nk}\cos{k\lambda} + b_{nk}\sin{k\lambda}\right) P_n^k (\cos\theta), $$ где $Y_n^k$ может быть вычислена по функции scipy.special.sph_harm.

Нормированные сферические функции

Как видно, в $Y_n^k (\theta, \lambda)$ присутствует нормирующий множитель $$ N_{nk} = \sqrt{\frac{(2n + 1)}{4\pi}\frac{(n - k)!}{(n + k)!}}, $$ поэтому сферические функции будут нормированными (normalized). Заметим, что они не будут полностью нормированными (fully normalized), ибо в этом случае нормирующий множитель должен быть $$ N_{nk} = \sqrt{(2n + 1)\frac{(n - k)!}{(n + k)!}}, $$ который легко получить из предыдущего выражения, умножив его на $\sqrt{4\pi}$, поэтому полностью нормированные сферические функции часто называют ещё $4\pi$ - normalized.

Фаза Кондона-Шортли

Под фазой Кондона-Шортли подразумевается множитель $(-1)^k$, который обычно вводится либо для присоединенных функций Лежандра, как это и сделано в SciPy для функции scipy.special.lpmv, либо для самих сферических функций. Применяется чаще всего в физике и сейсмологии. В геодезии этот множитель никогда не используется, поэтому от него необходимо избавиться, поделив сферическую функцию, вычисленную через scipy.special.sph_harm, на величину фазы Кондона-Шортли, то есть на $-1$ при нечетном $k$.

Пример вычисления сферической функции в SciPy

Вычислим в качестве примера, а заодно и проверим этим наши выкладки, сферическую функцию степени $n=5$ и порядка $k=3$. Импортируем всё необходимое

import numpy as np
from scipy.special import sph_harm, lpmv, factorial

и вычисляем комплексную сферическую функцию

n = 5
k = 3
# произвольные координаты в радианах
colat, lon = 0.3, 0.7
Y53 = sph_harm(k, n, lon, colat)
Y53

(0.03251642274478207-0.055598093017787176j)

Сделаем косвенное вычисление для проверки. Для начала получим значение присоединенной функции Лежандра $P_n^k (\cos\theta)$ с примененной по умолчанию фазой Кондона-Шортли. Сделать это можно через функцию scipy.special.lpmv так

x = np.cos(colat)
P53 = lpmv(k, n, x)
P53

-9.774568734291488

Теперь получим нормирующий множитель $N_{nk}$

Nnk = np.sqrt((2*n + 1) / (4 * np.pi) * factorial(n - k) / factorial(n + k))
Nnk

0.006589404174225528

Наконец, сравним встроенную реализацию функции scipy.special.sph_harm с полученной нами

Nnk * P53 * np.exp(k*lon*1j) == Y53

True

Результат одинаковый. Осталось посмотреть, правильно ли мы разобрались со связью вещественной и комплексной сферической функции. Для этого вычислим ненормированную сферическую функцию по привычному для геодезистов определению, уберём разве что фазу Кордона-Шортли, поскольку она автоматически добавляется в scipy.special.lpmv, а затем сравним результаты

# Произвольные вещественные и комплексный коэффициенты
ank, bnk = 0.5, 0.8
cnk = ank - bnk*1j

Y53_real = (ank * np.cos(k * lon) + bnk * np.sin(k * lon)) * P53 * (-1)**k

np.allclose((cnk * (Y53 / Nnk)).real * (-1)**k, Y53_real)

True

В конце стоит сказать об ограничениях. В настоящее время реализованный в SciPy алгоритм не позволяет вычислять присоединённые функции Лежандра высоких степений. Например, уже для $n = 86, k = 86$ имеем

print(lpmv(86, 86, x))
print(sph_harm(86, 86, lon, colat))

inf
(nan+nanj)

Это известная проблема почти всех библиотек для всех языков программирования. Для SciPy, например, уже больше года на GitHub висит не решённый Issue: Overflow in sph_harm for high degrees. Этот барьер в вычислениях накладывает сильные ограничения по использованию большинства доступных готовых решений для геодезических задач, где используемые степени в разложении явно больше. Но решение, конечно, имеется. Основная публикация здесь вот эта:

Holmes, S. A., and W. E. Featherstone, A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions, J. Geodesy, 76, 279- 299, doi:10.1007/s00190-002-0216-2, 2002.

Этот алгоритм реализован, например, в pyshtools, в которой можно свободно работать со сферическими функциями до 2800-й степени. Для простых задач, однако, например, для визуализации сферических функций низких степеней, возможностей SciPy вполне достаточно.