Conical: an extended module for computing a numerically satisfactory pair of solutions of the differential equation for conical functions

TL;DR

This paper introduces an extended computational module for conical functions, enabling accurate calculation of a new companion function crucial for solving boundary value problems in conical domains.

Contribution

The paper presents an extension to the CONICAL module, adding routines to compute a new numerically satisfactory conical function for improved problem-solving.

Findings

Provides a routine for computing ${R}^{m}_{-rac{1}{2}+i au}(x)$

Enables solving Dirichlet problems in conical domains

Enhances the computational tools for conical functions

Abstract

Conical functions appear in a large number of applications in physics and engineering. In this paper we describe an extension of our module CONICAL for the computation of conical functions. Specifically, the module includes now a routine for computing the function ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$, a real-valued numerically satisfactory companion of the function ${\rm P}^m_{-\tfrac12+i\tau}(x)$ for $x>1$. In this way, a natural basis for solving Dirichlet problems bounded by conical domains is provided.

Equations

$$(1-x^2)\frac{d^{2}w}{d{x}^2}-2x\frac{dw}{dx}-\left({\tau }^2+\frac{1}{4}+\frac{{\mu }^2}{1-x^2}\right)w=0.$$

$${\mathrm{P}}_{-\frac{1}{2}+i\tau }^{-m}(x)=\frac{\pi }{\cosh (\pi \tau ){\left|\Gamma (m+\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}+i\tau )\right|}^2}{\mathrm{P}}_{-\frac{1}{2}+i\tau }^m(x).$$

$$W\text{/}\left\{{\mathrm{P}}_{-\frac{1}{2}+i\tau }^m\text{/}\left(x\right),{\mathrm{R}}_{-\frac{1}{2}+i\tau }^m\text{/}\left(x\right)\right\}=\frac{\pi (e^{-\tau \pi }+\sinh \text{/}\left(\tau \pi \right))}{|\Gamma \text{/}\left(-m+\frac{1}{2}+i\tau \right){|}^2({\cosh \text{/}}^2\left(\tau \pi \right))(1-x^2)}.$$

$$\frac{d}{dx}{\mathrm{P}}_{-\frac{\text{lower}1.80835pt1}{\text{raise}1.50693pt2}+i\tau }^m(x)=-\frac{1}{\sqrt{x^2-1}}{\mathrm{P}}_{-\frac{\text{lower}1.80835pt1}{\text{raise}1.50693pt2}+i\tau }^{m+1}(x)+\frac{mx}{x^2-1}{\mathrm{P}}_{-\frac{\text{lower}1.80835pt1}{\text{raise}1.50693pt2}+i\tau }^m(x).$$

$$\begin{array}[]{lcl}{\rm R}_{-\frac{1}{2}+i\tau}^{\,0}(x)&=&\Re\left(\displaystyle\sum_{k=0}^{\infty}\frac{\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}-i\tau\right)_{k}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}+i\tau\right)_{k}}{(k!)^{2}}z^{k}\bigl{(}\psi(k+1)-\psi({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}+i\tau)-\ln w\bigr{)}\right),\\ {\rm R}_{-\frac{1}{2}+i\tau}^{\,1}(x)&=&\Re\left(\sqrt{x^{2}-1}\displaystyle\sum_{k=0}^{\infty}\frac{\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}-i\tau\right)_{k}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}+i\tau\right)_{k}}{(k!)^{2}}{{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}z^{k-1}\right.\ \left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}\left(-1+w^{2}\right)+\right.\\ &&\left.\left.k\bigl{(}\psi(k+1)-\psi({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}+i\tau)-\ln w\bigr{)}\right)\right),\end{array}$$

$$z=\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}(1-x),w=\sqrt{\frac{x-1}{x+1}},$$

$$\psi (\alpha )\sim \ln \alpha -\frac{1}{2\alpha }-\frac{1}{12{\alpha }^2}+\frac{1}{120{\alpha }^4}-\frac{1}{252{\alpha }^6}+...,$$

$${\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m+1}(x)-\frac{2mx}{\sqrt{x^2-1}}{\mathrm{R}}_{-\frac{1}{2}+i\tau }^m(x)+\left((m-\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}{)}^2+{\tau }^2\right){\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m-1}(x)=0$$

$$\begin{array}[]{@{}r@{\;}c@{\;}l@{}}{\rm R}_{-\frac{1}{2}+i\tau}^{\mu}(x)&\sim&\displaystyle\sqrt{\pi/2}\,\alpha^{\mu+\frac{1}{2}}\left(x^{2}-1\right)^{-\frac{1}{4}}\ \times\\[8.0pt] &&\displaystyle\sum_{k=0}^{\infty}\,f_{k}\left(\Re\Phi_{k}\cos\phi+\Im\Phi_{k}\sin\phi\right),\end{array}$$

$$z=\frac{1}{2\sqrt{x^2-1}\left(x+\sqrt{x^2-1}\right)}.$$

$${\Phi }_k={\left(\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}-\mu \right)}_k{\omega }^{2\mu -k}U\left(\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}+\mu ,1+2\mu -k,\alpha \omega \right),\omega =i\tau .$$

$$H_{\mu }^{(2)}(z)=J_{\mu }(z)-i{Y}_{\mu }(z).$$

$$\begin{array}[]{lcl}\Phi_{0}&=&-{{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}i\sqrt{\pi}(\tau/\alpha)^{\mu}e^{\frac{1}{2}i\alpha\tau}H_{\mu}^{(2)}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}\alpha\tau\right),\\ \Phi_{1}&=&\displaystyle{{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 4$}}}\alpha\sqrt{\pi}(\tau/\alpha)^{\mu}e^{\frac{1}{2}i\alpha\tau}\left(iH_{\mu}^{(2)}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}\alpha\tau\right)+H_{\mu-1}^{(2)}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}\alpha\tau\right)\right).\end{array}$$

$$\omega {\Phi }_{n+1}=(n-2\mu -\alpha \omega ){\Phi }_n+\alpha (n-\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}-\mu ){\Phi }_{n-1}.$$

$$\begin{array}[]{lcl}\Re\Phi_{n+1}&=&\frac{\displaystyle{n-2\mu}}{\displaystyle{\tau}}\Im\Phi_{n}-\alpha\Re\Phi_{n}+\frac{\displaystyle{\alpha}}{\displaystyle{\tau}}(n-{{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}-\mu)\Im\Phi_{n-1},\\ \Im\Phi_{n+1}&=&-\frac{\displaystyle{n-2\mu}}{\displaystyle{\tau}}\Re\Phi_{n}-\alpha\Im\Phi_{n}-\frac{\displaystyle{\alpha}}{\displaystyle{\tau}}(n-{{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 2.1097pt\hbox{$\scriptstyle 2$}}}-\mu)\Re\Phi_{n-1}.\end{array}$$

$$\begin{array}[]{lcl}f_{0}&=&1,\\[6.0pt] f_{1}&=&\frac{\displaystyle{b}}{\displaystyle{2d}}\left(2dz+d-2z\right),\\[6.0pt] f_{2}&=&\frac{\displaystyle{b}}{\displaystyle{24d^{2}}}\left(12z^{2}+12bz^{2}+d^{2}-12d^{2}z-12d^{2}z^{2}-24bdz^{2}\right.\\[6.0pt] &&\left.\,\,\,\,\,\,\,\,\,\,+12bd^{2}z+12bd^{2}z^{2}+3bd^{2}-12bdz\right),\end{array}$$

$$\begin{array}[]{@{}r@{\;}c@{\;}l@{}}{\rm R}_{-\frac{1}{2}+i\tau}^{\mu}(x)&=&\Re\left(\displaystyle\sqrt{\pi/2}\,\left(x^{2}-1\right)^{-\frac{1}{4}}G(\mu,\tau)\,e^{-i\phi}\right.\ \times\\ &&\left.\displaystyle\sum_{k=0}^{\infty}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}+\mu\right)_{k}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}-\mu\right)_{k}\frac{u_{k}(\tau)+iv_{k}(\tau)}{w_{k}(\tau)}\frac{\displaystyle{(-z)^{k}}}{\displaystyle{k!}}\right),\end{array}$$

$$G(\mu ,\tau )=\frac{\Gamma \left(\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}+\mu +i\tau \right)}{\Gamma \left(1+i\tau \right)},$$

$$\frac{u_k(\tau )+i{v}_k(\tau )}{w_k(\tau )}=\frac{1}{{\left(1+i\tau \right)}_k},k=0,1,2,\dots .$$

$$\begin{array}[]{@{}r@{\;}c@{\;}l@{}}u_{k+1}(\tau)&=&(k+1)u_{k}(\tau)+\tau v_{k}(\tau),\\[8.0pt] v_{k+1}(\tau)&=&(k+1)v_{k}(\tau)-\tau u_{k}(\tau),\\[8.0pt] w_{k+1}(\tau)&=&\left((k+1)^{2}+\tau^{2}\right)w_{k}(\tau),\end{array}$$

$$G(\mu ,\tau )=H(\mu ,\tau )e^{i\rho (\mu ,\tau )},u_k(\tau )+i{v}_k(\tau )=r_k(\tau )e^{i{\sigma }_k(\tau )},$$

$$\begin{array}[]{lcl}{\rm R}_{-\frac{1}{2}+i\tau}^{\mu}(x)&=&\displaystyle\sqrt{\pi/2}\,H(\mu,\tau)\,\left(x^{2}-1\right)^{-\frac{1}{4}}\,\ \times\\ &&\displaystyle\sum_{k=0}^{\infty}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}+\mu\right)_{k}\left({{\lower 2.58334pt\hbox{$\scriptstyle 1$}}\over{\raise 3.01385pt\hbox{$\scriptstyle 2$}}}-\mu\right)_{k}\frac{\displaystyle{r_{k}(\tau)}}{\displaystyle{w_{k}(\tau)}}\frac{\displaystyle{(-z)^{k}}}{\displaystyle{k!}}\cos(\psi_{k}),\end{array}$$

$${\psi }_k=\tau \log \left(x+\sqrt{x^2-1}\right)-\rho (\mu ,\tau )-{\sigma }_k(\tau ).$$

$$\frac{\left(2mx/\sqrt{x^2-1}\right){\mathrm{R}}_{-\frac{1}{2}+i\tau }^m(x)-\left((m-\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}{)}^2+{\tau }^2\right){\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m-1}(x)}{{\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m+1}(x)}=1.$$

$$\frac{{\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m+1}(x)+\left((m-\frac{\text{lower}2.58334pt1}{\text{raise}3.01385pt2}{)}^2+{\tau }^2\right){\mathrm{R}}_{-\frac{1}{2}+i\tau }^{m-1}(x)}{\left(2mx/\sqrt{x^2-1}\right){\mathrm{R}}_{-\frac{1}{2}+i\tau }^m(x)}=1.$$

$${\mathrm{P}}_{-\frac{1}{2}+i\tau }^1(x){\mathrm{R}}_{-\frac{1}{2}+i\tau }^0(x)-{\mathrm{P}}_{-\frac{1}{2}+i\tau }^0(x){\mathrm{R}}_{-\frac{1}{2}+i\tau }^1(x)=\frac{e^{-\pi \tau }+\sinh (\pi \tau )}{\cosh (\pi \tau )\sqrt{x^2-1}},$$

Full text

The module also improves the performance of our previous algorithm for the conical function {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) and it includes now the computation of the first order derivative of the function. This is also considered for the function ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$ in the extended algorithm.

††journal: Computer Physics Communications

1 Introduction

Conical or Mehler functions are involved in a large number of applications in different areas of physics. In particular, these functions appear when solving the Laplace equation in spherical coordinates for two intersecting cones [2] or for regions bounded by two intersecting spheres, or by one or two confocal hyperboloids of revolution when using toroidal coordinates.

An extended version of our module Conical [1] for the computation of conical functions is presented in this paper. The new module includes now a routine for computing the function ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$, a real-valued numerically satisfactory companion of the function {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) for $x>1$. The module also improves our previous algorithm for the conical function {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) by considering now more cofficients in some of the asymptotic expansions used for computing the function in the region $x>1$. The computation of the first order derivatives of {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) and ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$ is also included in the new module.

2 Theoretical background

Conical functions are solutions of the second order differential equation

[TABLE]

We will restrict to integer positive values of the parameter $\mu$ ($\mu=m\in{\mathbb{Z}}^{+}$).

When $-1<x<1$, a real-valued satisfactory pair of solutions of (1) is ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm P}^{-m}_{-\frac{1}{2}+i\tau}(x)$. Both functions can be computed using our algorithm for ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ implemented in conicp; for computing ${\rm P}^{-m}_{-\frac{1}{2}+i\tau}(x)$ the following relation can be used

[TABLE]

When $x>1$, a real-valued satisfactory pair of solutions of (1) is ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)\equiv\Re\left\{e^{-i\pi m}{{\rm Q}}^{m}_{-\frac{1}{2}+i\tau}(x)\right\}$ (the function ${{\rm Q}}^{m}_{-\frac{1}{2}+i\tau}(x)$ is complex-valued).

The Wronskian relation for ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$, useful for testing, is given by

[TABLE]

The algorithm for computing the conical function ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ was described in [1]. In the new module Conical we also compute the first order derivatives for {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) and ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$: the first order derivative of {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) can be obtained using the relation

[TABLE]

The first order derivative of ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$ satisfies the same relation (4) with {\rm R}^{m+1}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) and {\rm R}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x).

A preliminary algorithm for computing the function ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ was presented in [3] although the final algorithm in finite precision arithmetic implemented in the routine conicr presents some differences with respect to the first algorithm. We have also changed the notation of the function with respect to the one used in that reference; we are using now ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ instead of $\widetilde{Q}^{m}_{-\frac{1}{2}+i\tau}(x)$ for simplicity.

Next we summarize the methods for computing the conical function ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$.

2.1 Computation of ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ for values of $x$ close to $1$

2.1.1 Small or moderate values of $\tau$

To compute ${\rm R}^{0}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{1}_{-\frac{1}{2}+i\tau}(x)$ we will use the expansions

[TABLE]

where $z$ and $w$ are given by

[TABLE]

and $\psi(\alpha)=\Gamma^{\prime}(\alpha)/\Gamma(\alpha)$.

The only complex quantity in (5) is the function \psi({{\lower 2.58334pt\hbox{\scriptstyle 1}}\over{\raise 2.1097pt\hbox{\scriptstyle 2}}}+i\tau). For computing this function we use an algorithm which computes separately the real and imaginary parts of \psi({{\lower 2.58334pt\hbox{\scriptstyle 1}}\over{\raise 2.1097pt\hbox{\scriptstyle 2}}}+i\tau) avoiding, in this way, the use of complex arithmetics when computing (5). The algorithm is based on the use of the asymptotic expansion

[TABLE]

valid for $\alpha\to\infty$ in $|{\rm ph}\,\alpha|<\pi$. We use this expansion if $|\alpha|\geq 12$ with $8$ terms of the series (or less), and we use the recurrencence relation $\psi(\alpha)=\psi(\alpha+1)-1/\alpha$ for smaller values of $|\alpha|$.

When $m\geq 2$ we will use the recursion relation

[TABLE]

in the direction of increasing $m$.

2.1.2 Large values of $\tau$

A representation in terms of Kummer $U-$functions is used in this case:

[TABLE]

where $\phi=\tau\log(x+\sqrt{x^{2}-1})$, $\alpha=\ln\left(\frac{\displaystyle{z+1}}{\displaystyle{z}}\right)$ and $z$ is given by

[TABLE]

The functions $\Phi_{k}$ are given in terms of Kummer $U-$functions as follows

[TABLE]

The functions can be also written in terms of the Hankel functions $H_{\mu}^{(2)}(z)$. This representation makes simple separating the real and imaginary parts of $\Phi_{k}$ by using

[TABLE]

For the computation of the Bessel functions $J_{\mu}(z)$, $Y_{\mu}(z)$ we use an algorithm which combines the use of series expansions, Debye’s asymptotic expansions, asymptotic expansions for large $z$, Airy-type asymptotic expansions and three-term recurrence relations. This algorithm is implemented in the module BesselJY and it is also included in the software package.

The functions $\Phi_{0}$, $\Phi_{1}$ are given by

[TABLE]

For computing $\Phi_{n}$ for $n=2,...$ we can use a recurrence relation for the Kummer $U-$functions which gives

[TABLE]

From this, the following recurrence relations for both the real and imaginary parts of $\Phi_{n+1}$ can be obtained:

[TABLE]

The first coefficients $f_{k}$ in (9) are given by

[TABLE]

where $b=-\mu-\frac{1}{2}$ and $d=z\alpha$.

2.2 Computation of ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ for moderate or large values of $x$

In this case we use the expansion

[TABLE]

where $z$ is given in eq.(10), $\phi=\tau\log\left(x+\sqrt{x^{2}-1}\right)$,

[TABLE]

and

[TABLE]

We can compute $u_{k}(\tau)$, $v_{k}(\tau)$ and $w_{k}(\tau)$ from the recurrence relations

[TABLE]

with $u_{0}(\tau)=1$, $v_{0}(\tau)=0$, $w_{0}(\tau)=1$.

The real part of (17) can be obtain by writing

[TABLE]

which gives

[TABLE]

where

[TABLE]

The computation of the ratio of two gamma functions in (18) is made using an algorithm for computing the gamma function for complex values of the argument. The algorithm adapts for complex arguments the scheme used for real values described in [4].

3 Overview of the software structure

The Fortran 90 package includes the main module Conical, which includes the routines conicp, conicr and conicpr.

In the module Conical, the auxiliary module Someconstants is used. This is a module for the computation of the main constants used in the different routines. The module BesselJY (for the computation of Bessel functions) and AiryFunction (for the computation of Airy functions) are used. The routines included in auxil.f90 are also used in the module Conical.

4 Description of the individual software components

The Fortran 90 module Conical includes the public routine conicp, which computes the conical functions ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$ for $x>-1$, $m\geq 0$ and $\tau>0$; the routine conicr, which computes the function ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$, for $x>1$ , $m\geq 0$ and $\tau>0$ and the routine conicpr, which computes both functions ${\rm P}^{m}_{-\frac{1}{2}+i\tau}(x)$, ${{\rm R}}^{m}_{-\frac{1}{2}+i\tau}(x)$ and their first order derivatives for $x>1$ , $m\geq 0$ and $\tau>0$. The calling sequences of these routines are

CALL conicp(x,mu,tau,pm,ierr)

where the input data are: $x$, $mu$ and $tau$ (arguments of the functions). The outputs are the error flag $ierr$ and the function value $pm$. The possible values of the error flag are: $ierr=0$, successful computation; $ierr=1$, computation failed due to overflow/underflow; $ierr=2$, arguments out of range.

CALL conicr(x,mu,tau,rm,ierr)

where the input data are: $x$, $mu$ and $tau$ (arguments of the functions). The outputs are the error flag $ierr$ and the function value $rm$. The possible values of the error flag are: $ierr=0$, successful computation; $ierr=1$, computation failed due to overflow/underflow; $ierr=2$, arguments out of range.

CALL conicpr(x,mu,tau,pm,pmd,rm,rmd,ierr)

where the input data are: $x$, $mu$ and $tau$ (arguments of the functions). The outputs are the error flag $ierr$, the function values $pm$, $rm$ and the first order derivatives $pmd$, $rmd$. The possible values of the error flag are: $ierr=0$, successful computation; $ierr=1$, computation failed.

5 Testing the algorithm

For testing the accuracy of the expansions used to compute the conical function ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$, we have first checked (8) written in the form

[TABLE]

This check fails close to the zeros of ${\rm R}^{m+1}_{-\frac{1}{2}+i\tau}(x)$; in this case, we can consider the alternative test

[TABLE]

Because the zeros of ${\rm R}^{m+1}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ are interlaced, both tests will not fail simultaneously. We can therefore take the minimum of both errors. We have considered these tests for the expansions described in sections 2.1.2 ($x$ small) and 2.2 ($x$ large). Figure 1 shows the points where the minimum value of the error of the tests (24) and (25) when using (9) is greater than $10^{-12}$. In the algorithm, we have fixed to $N=7$ the number of terms used in the expansion. Random points have been generated in the domain $(x,\,\tau)\in(1.001,\,1.05)\times(15,\,100)$. As can be seen, for $m=1$ (upper figure) the use of the expansion (9) allows to compute the function values with an accuracy better than $10^{-12}$ for values of $\tau$ greater than $20$ when $x$ is close to $1$. The accuracy of the expansion worsens as $m$ increases, as can be seen also in Figure 1 (lower figure) where the same test is considered for $m=5$. Therefore, one has to use an alternative method of computation for moderate/large values of $m$ as, for example, the use of the recurrence relation (8) starting from ${\rm R}^{0}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{1}_{-\frac{1}{2}+i\tau}(x)$.

Figure 2 shows the same tests (24) and (25) for the expansion (17) and for $\mu\equiv m=1$. The domain where the random points have been generated is now $(x,\,\tau)\in(1.2,\,100)\times(0,\,100)$. As can be seen in the figure, there is some loss of accuracy when $\tau$ is large and $x$ is moderate/large. In any case, we have checked that the accuracy was always better than $5\,10^{-12}$.

For testing the expansions for ${\rm R}^{0}_{-\frac{1}{2}+i\tau}(x)$ and ${\rm R}^{1}_{-\frac{1}{2}+i\tau}(x)$ of section 2.1.1, we have used the Wronskian relation given in (3). In this case, we have

[TABLE]

Figure 3 shows the points where the value of the error in the Wronskian relation (26) is greater than $10^{-12}$.

The accuracy of the final algorithm for ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ has been tested by computing the Wronskian relation given in (3) for a very large number of random points in the parameter domain $(x,\,m\,\tau)\in(1.001,\,100)\times(0,\,100)\times(0,\,100)$. The algorithm for the conical function {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x) was improved by considering more cofficients in some of the asymptotic expansions used for computing the function in the region $x>1$. We have checked that the accuracy of the Wronskian test (3) is close to $10^{-12}$ in the whole parameter domain and better than $10^{-13}$ for a large fraction of the tested parameter values.

6 Test run description

The Fortran 90 test program testcon.f90 includes the computation of 20 values of the functions {\rm P}^{m}_{-{{\lower 1.80835pt\hbox{\scriptstyle 1}}\over{\raise 1.50693pt\hbox{\scriptstyle 2}}}+i\tau}(x), ${\rm R}^{m}_{-\frac{1}{2}+i\tau}(x)$ and their first order derivatives and their comparison with the corresponding pre-computed results.

7 Acknowledgements

A.G. acknowledges the Fulbright/MEC Program for support during her stay at SDSU. J.S. acknowledges the Salvador de Madariaga Program for support during his stay at SDSU. The authors acknowledge financial support from Ministerio de Ciencia e Innovación, project MTM2015-67142-P. NMT thanks CWI, Amsterdam, for scientific support.