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

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.


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 R m − 1 2 +iτ (x), a real-valued numerically satisfactory companion of the function P m − 1 2 +iτ (x) for x > 1.The module also improves our previous algorithm for the conical function P m − 1 2 +iτ (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 P m − 1 2 +iτ (x) and R m − 1 2 +iτ (x) is also included in the new module.

Theoretical background
Conical functions are solutions of the second order differential equation We will restrict to integer positive values of the parameter µ (µ = m ∈ Z + ).
When −1 < x < 1, a real-valued satisfactory pair of solutions of (1) is P m − 1 2 +iτ (x) and P −m − 1 2 +iτ (x).Both functions can be computed using our algorithm for P m − 1 2 +iτ (x) implemented in conicp; for computing P −m − 1 2 +iτ (x) the following relation can be used When x > 1, a real-valued satisfactory pair of solutions of ( 1) is P m The Wronskian relation for P m − 1 2 +iτ (x) and R m − 1 2 +iτ (x), useful for testing, is given by .
The algorithm for computing the conical function P m − 1 2 +iτ (x) was described in [1].In the new module Conical we also compute the first order derivatives for P m − 1 2 +iτ (x) and R m − 1 2 +iτ (x): the first order derivative of P m − 1 2 +iτ (x) can be obtained using the relation The first order derivative of R m A preliminary algorithm for computing the function R m − 1 2 +iτ (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 R m where z and w are given by and ψ(α) = Γ ′ (α)/Γ(α).
The only complex quantity in (5) is the function ψ( 1 2 + iτ ).For computing this function we use an algorithm which computes separately the real and imaginary parts of ψ( 12 + iτ ) avoiding, in this way, the use of complex arithmetics when computing (5).The algorithm is based on the use of the asymptotic expansion valid for α → ∞ in |ph α| < π.We use this expansion if |α| ≥ 12 with 8 terms of the series (or less), and we use the recurrencence relation When m ≥ 2 we will use the recursion relation in the direction of increasing m.

Large values of τ
A representation in terms of Kummer U −functions is used in this case: where The functions Φ k are given in terms of Kummer U −functions as follows The functions can be also written in terms of the Hankel functions H (2) µ (z).This representation makes simple separating the real and imaginary parts of Φ k by using For the computation of the Bessel functions J µ (z), Y µ (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 Φ 0 , Φ 1 are given by For computing Φ n for n = 2, ... we can use a recurrence relation for the Kummer U −functions which gives From this, the following recurrence relations for both the real and imaginary parts of Φ n+1 can be obtained: The first coefficients f k in (9) are given by where b = −µ − 1 2 and d = zα.

Computation of R m
− 1 2 +iτ (x) for moderate or large values of x In this case we use the expansion where z is given in eq.( 10), and We can compute u k (τ ), v k (τ ) and w k (τ ) from the recurrence relations with u 0 (τ ) = 1, v 0 (τ ) = 0, w 0 (τ ) = 1.The real part of (17) can be obtain by writing which gives where 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].

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.f90are also used in the module Conical.

Description of the individual software components
The Fortran 90 module Conical includes the public routine conicp, which computes the conical functions P m 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.

Testing the algorithm
For testing the accuracy of the expansions used to compute the conical function R m − 1 2 +iτ (x), we have first checked (8) written in the form This check fails close to the zeros of R m+1 − 1 2 +iτ (x); in this case, we can consider the alternative test Because the zeros of R m+1 − 1 2 +iτ (x) and R m − 1 2 +iτ (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, τ ) ∈ (1.001, 1.05) × (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 τ 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 R 0    Figure 2 shows the same tests ( 24) and (25) for the expansion (17) and for µ ≡ m = 1.The domain where the random points have been generated is now (x, τ ) ∈ (1.2, 100) × (0, 100).As can be seen in the figure, there is some loss of accuracy when τ 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 R 0 − 1 2 +iτ (x) and R 1 − 1 2 +iτ (x) of section 2.1.1,we have used the Wronskian relation given in (3).In this case, we have 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 R m − 1 2 +iτ (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 τ ) ∈ (1.001, 100) × (0, 100) × (0, 100).The algorithm for the conical function P m − 1 2 +iτ (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.

Figure 1 :
Figure 1: Test of the performance of the expansion (9).The points where the value of the error when testing the recurrence relation (8) is greater than 10 −12 are plotted.

Figure 2 :
Figure 2: Test of the performance of the expansion (17) for µ ≡ m = 1.The points where the value of the error in the recurrence relation (8) is greater than 10 −12 are plotted.

Figure 3 :
Figure 3: Test of the performance of the expansions given in (5).The points where the value of the error in the Wronskian relation (26) is greater than 10 −12 are plotted.