implement spherical harmonics

acb_hypgeom.h

@@ -169,6 +169,7 @@ void acb_hypgeom_laguerre_l(acb_t res, const acb_t n, const acb_t m, const acb_t
 void acb_hypgeom_hermite_h(acb_t res, const acb_t n, const acb_t z, long prec);
 void acb_hypgeom_chebyshev_t(acb_t res, const acb_t n, const acb_t z, long prec);
 void acb_hypgeom_chebyshev_u(acb_t res, const acb_t n, const acb_t z, long prec);
+void acb_hypgeom_spherical_y(acb_t res, long n, long m, const acb_t theta, const acb_t phi, long prec);
 
 #ifdef __cplusplus
 }

acb_hypgeom/spherical_y.c

@@ -0,0 +1,93 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_spherical_y(acb_t res, long n, long m,
+    const acb_t theta, const acb_t phi, long prec)
+{
+    acb_t t, u;
+
+    if (n < 0)
+    {
+        if (m <= n)
+        {
+            acb_zero(res);
+            return;
+        }
+
+        n = -1-n;
+    }
+
+    if (m > n || m < -n)
+    {
+        acb_zero(res);
+        return;
+    }
+
+    if (n > LONG_MAX / 4)
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    acb_init(t);
+    acb_init(u);
+
+    acb_sin_cos(t, u, theta, prec);
+
+    /* P_n^m(cos(theta)) */
+    acb_hypgeom_legendre_p_uiui_rec(u, n, FLINT_ABS(m), u, prec);
+    acb_pow_ui(t, t, FLINT_ABS(m), prec);
+    acb_mul(t, t, u, prec);
+
+    /* exp(i m phi) */
+    acb_mul_onei(u, phi);
+    acb_mul_si(u, u, m, prec);
+    acb_exp(u, u, prec);
+    if (m < 0 && m % 2)
+        acb_neg(u, u);
+    acb_mul(t, t, u, prec);
+
+    /* sqrt((2n+1)/(4pi) (n-m)!/(n+m)!) */
+    arb_fac_ui(acb_realref(u), n - FLINT_ABS(m), prec);
+    arb_fac_ui(acb_imagref(u), n + FLINT_ABS(m), prec);
+    arb_mul_ui(acb_realref(u), acb_realref(u), 2 * n + 1, prec);
+    arb_div(acb_realref(u), acb_realref(u), acb_imagref(u), prec);
+
+    arb_const_pi(acb_imagref(u), prec);
+    arb_div(acb_realref(u), acb_realref(u), acb_imagref(u), prec);
+    arb_mul_2exp_si(acb_realref(u), acb_realref(u), -2);
+    arb_sqrt(acb_realref(u), acb_realref(u), prec);
+
+    acb_mul_arb(t, t, acb_realref(u), prec);
+
+    acb_set(res, t);
+
+    acb_clear(t);
+    acb_clear(u);
+}
+

acb_hypgeom/test/t-spherical_y.c

@@ -0,0 +1,167 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+/* generated with mpmath */
+static const double testdata[] = {
+    0.0, 0.0,
+    -0.02701626593453869896, 0.06085851132589568114,
+    0.092391418427377575111, 0.29443833672801130076,
+    0.65131374739365937388, 0.37063665572827331764,
+    0.83847032733386007862, -0.28730577566514336534,
+    -0.14750357092078405673, -0.30269451746476217073,
+    -0.048646822931479694311, 0.038926542608778397121,
+    0.0048197328091099892156, 0.003641138129582944611,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.027574808444335508972, 0.11025410152651193448,
+    0.35400136349600893569, 0.24420566191536183849,
+    0.67425312681608106523, -0.11728998324910652574,
+    -0.069322941709900597259, -0.18037674817459565026,
+    -0.018729775699827662658, 0.013254956213452933465,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.14771800806690767222, 0.11532328561497280466,
+    0.50057392828327004059, -0.029559886614345590592,
+    -0.02552211332726607337, -0.080244885187214384683,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.28209479177387814347, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.28209479177387814347, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.14771800806690767222, 0.11532328561497280466,
+    0.50057392828327004059, -0.029559886614345590592,
+    -0.02552211332726607337, -0.080244885187214384683,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.027574808444335508972, 0.11025410152651193448,
+    0.35400136349600893569, 0.24420566191536183849,
+    0.67425312681608106523, -0.11728998324910652574,
+    -0.069322941709900597259, -0.18037674817459565026,
+    -0.018729775699827662658, 0.013254956213452933465,
+    0.0, 0.0,
+    0.0, 0.0,
+    0.0, 0.0,
+    -0.02701626593453869896, 0.06085851132589568114,
+    0.092391418427377575111, 0.29443833672801130076,
+    0.65131374739365937388, 0.37063665572827331764,
+    0.83847032733386007862, -0.28730577566514336534,
+    -0.14750357092078405673, -0.30269451746476217073,
+    -0.048646822931479694311, 0.038926542608778397121,
+    0.0048197328091099892156, 0.003641138129582944611,
+    0.0, 0.0,
+    -0.033798002071615353821, 0.018033964330683617275,
+    -0.071988993423855434033, 0.19195223773954952294,
+    0.23381678888295549361, 0.58288246446413288682,
+    1.0666048448920927459, 0.45816450098677839043,
+    0.99160821344691798668, -0.56670953266315770842,
+    -0.27930640165163989333, -0.44051756827083040505,
+    -0.092578419347396006256, 0.086641542063743427664,
+    0.015474322098595795308, 0.010316287906840124695,
+    0.00051935757617903067671, -0.0014726339944978705874,
+};
+
+int main()
+{
+    printf("spherical_y....");
+    fflush(stdout);
+
+    {
+        long i, n, m;
+        acb_t z, w, x, y;
+
+        acb_init(z);
+        acb_init(w);
+        acb_init(x);
+        acb_init(y);
+
+        i = 0;
+
+        arb_set_str(acb_realref(x), "0.2", 64);
+        arb_set_str(acb_imagref(x), "0.3", 64);
+        arb_set_str(acb_realref(y), "0.3", 64);
+        arb_set_str(acb_imagref(y), "0.4", 64);
+
+        for (n = -4; n <= 4; n++)
+        {
+            for (m = -4; m <= 4; m++)
+            {
+                acb_hypgeom_spherical_y(z, n, m, x, y, 64);
+
+                acb_set_d_d(w, testdata[2 * i], testdata[2 * i + 1]);
+                mag_set_d(arb_radref(acb_realref(w)), 1e-13);
+                mag_set_d(arb_radref(acb_imagref(w)), 1e-13);
+
+                if (!acb_overlaps(z, w))
+                {
+                    printf("FAIL: value\n\n");
+                    printf("n = %ld, m = %ld\n", n, m);
+                    printf("z = "); acb_printd(z, 20); printf("\n\n");
+                    printf("w = "); acb_printd(w, 20); printf("\n\n");
+                    abort();
+                }
+
+                i++;
+            }
+        }
+
+        acb_clear(z);
+        acb_clear(w);
+        acb_clear(x);
+        acb_clear(y);
+    }
+
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_hypgeom.rst

@@ -787,3 +787,17 @@ Orthogonal polynomials and functions
     For nonnegative integer *n* and *m*, uses recurrence relations to evaluate
     `(1-z^2)^{-m/2} P_n^m(z)` which is a polynomial in *z*.
 
+.. function:: void acb_hypgeom_spherical_y(acb_t res, long n, long m, const acb_t theta, const acb_t phi, long prec)
+
+    Computes the spherical harmonic of degree *n*, order *m*,
+    latitude angle *theta*, and longitude angle *phi*, normalized
+    such that
+
+    .. math ::
+
+        Y_n^m(\theta, \phi) = \sqrt{\frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!}} e^{im\phi} P_n^m(\cos(\theta)).
+
+    The definition is extended to negative *m* and *n* by symmetry.
+    This function is a polynomial in `\cos(\theta)` and `\sin(\theta)`.
+    We evaluate it using :func:`acb_hypgeom_legendre_p_uiui_rec`.
+