yrecursion.f

cc Copyright (C) 2009-2012: Leslie Greengard and Zydrunas Gimbutas
cc Contact: greengard@cims.nyu.edu
cc
cc This software is being released under a modified FreeBSD license
cc (see COPYING in home directory).
cc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc License reproduced below:
cc
cc Copyright (c) 2009-2012, Leslie Greengard, Zydrunas Gimbutas
cc All rights reserved.
cc
cc Redistribution and use in source and binary forms, with or without
cc modification, are permitted provided that the following conditions are met:
cc
cc 1. Redistributions of source code must retain the above copyright notice, this
cc    list of conditions and the following disclaimer.
cc
cc 2. Redistributions in binary form must reproduce the above copyright
cc    notice, this list of conditions and the following disclaimer in the
cc    documentation and/or other materials provided with the distribution.
cc
cc 3. Neither the name of the copyright holder nor the names of its
cc    contributors may be used to endorse or promote products derived
cc    from this software without specific prior written permission.
cc
cc THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
cc ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
cc WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
cc DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
cc ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
cc (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
cc LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
cc ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
cc (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
cc SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    $Date$
c    $Revision$
c
c
c  This file contains a suite of evaluation codes for associated
c  Legendre functions with various scalings, arguent types, etc.
c  Following is a brief description of the user-callable subroutines.
c  (FORTRAN 77 VERSION).
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Simple routines for real argument.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndr - evaluate normalized Legendre functions
c
c       ylgndr2 - evaluate normalized Legendre functions
c                 and their derivatives
c
c       ylgndr2s - evaluate normalized Legendre functions and their
c            derivatives with scaling.
c            For m>0, the Ynm(x) values are scaled by 1/sqrt(1-x^2),
c            For m>0, the Ynm(x) derivatives with respect to x are
c            scaled by sqrt(1-x^2).
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Simple routines for real argument (not scaled by sqrt(2*n+1))
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndru - evaluate normalized Legendre functions
c
c       ylgndru2 - evaluate normalized Legendre functions
c                 and their derivatives
c
c       ylgndru2s - evaluate normalized Legendre functions and their
c            derivatives with scaling.
c            For m>0, the Ynm(x) values are scaled by 1/sqrt(1-x^2),
c            For m>0, the Ynm(x) derivatives with respect to x are
c            scaled by sqrt(1-x^2).
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Simple routines for real argument
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndr2sm - evaluate normalized Legendre functions and their
c            derivatives with scaling.
c            For m>0, the Ynm(x) values are scaled by 1/sqrt(1-x^2)**m,
c            For m>0, the Ynm(x) derivatives with respect to x are
c            scaled by 1/sqrt(1-x^2)**(m-2).
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Accelerated codes with precompution for real argument.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndrini - precomputation routine for recursion coefficients.
c                   Used by "fast" routines ylgndrf, ylgndr2f, ylgndr2sf.
c
c       ylgndrf - faster evaluation of normalized Legendre functions.
c                 Requires prior call to ylgndrini.
c
c       ylgndr2f - faster evaluation of normalized Legendre functions
c                  and their derivatives.
c                  Requires prior call to ylgndrini.
c
c       ylgndr2sf - faster evaluation of  normalized Legendre functions
c                   and their derivatives (with scaling as in ylgndr2s).
c                   Requires prior call to ylgndrini.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Simple codes for complex argument.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       zylgndr - evaluate normalized Legendre functions
c                 of a complex argument
c
c       zylgndr2 - evaluate normalized Legendre functions
c                  of a complex argument and their derivatives
c
c       zylgndr2s - evaluate normalized Legendre functions
c                  of a complex argument and their derivatives
c                  (with scaling as in ylgndr2s)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Accelerated codes with precompution for complex argument.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       zylgndrf - faster evaluation of normalized Legendre functions
c                  of a complex argument.
c                  Requires prior call to ylgndrini.
c
c       zylgndr2f - faster evaluation of normalized Legendre functions
c                   of a complex argument and their derivatives.
c                   Requires prior call to ylgndrini.
c
c       zylgndr2sf - faster evaluation of normalized Legendre functions
c                    of a complex argument and their derivatives
c                   (with scaling as in ylgndr2s).
c                   Requires prior call to ylgndrini.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Miscellaneous, special purpose codes.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       zylgndrbr - evaluate normalized Legendre functions of a
c            complex argument with modified branch cut.
c
c       zylgndrsc - evaluate normalized Legendre functions of a complex
c            argument with scaling to prevent overflow for large values
c            of z.
c
c       ylgndrpm, ylgndrpm_opt - Given Y_nm(x), return Y_nm(-x)
c
c       ylgndr2pm. ylgndrpm_opt - Given Y_nm(x) and Y'_nm(x),
c                                 return Y_nm(-x) and Y'_nm(-x)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Truncated evaluation routines: computation carried out for
c       orders up to specified parameter rather than max degree.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndr2s_trunc - evaluate normalized Legendre functions and their
c            derivatives (with scaling as in ylgndr2s),
c            BUT recursion carried out only to m = m2 rather than nmax.
c
c       ylgndrf_trunc - faster evaluation of normalized Legendre functions.
c            Same as ylgndrf, BUT recursion carried out only to m = m2
c            rather than nmax.
c            Requires prior call to ylgndrini.
c
c       ylgndr2f_trunc - faster evaluation of  normalized Legendre
c            functions and their derivatives. Same as ylgndr2f,
c            BUT recursion carried out only to m = m2 rather than nmax.
c
c       ylgndr2sf_trunc - faster evaluation of normalized Legendre
c            functions and their derivatives (with scaling).
c            Same as ylgndr2sf, BUT recursion carried out only to
c            m = m2 rather than nmax.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Accelerated codes with precompution for real argument:
c       precomputation is done only once for orders up to nmax,
c       the routines are able to revert to simple routines for
c       higher orders
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       ylgndrfwini - precomputation routine for recursion coefficients.
c                   Used by "fast" routines ylgndrfw, ylgndr2fw, ylgndr2sfw.
c
c       ylgndrfw - faster evaluation of normalized Legendre functions.
c                 Requires prior call to ylgndrfwini.
c
c       ylgndr2fw - faster evaluation of normalized Legendre functions
c                  and their derivatives.
c                  Requires prior call to ylgndrfwini.
c
c       ylgndr2sfw - faster evaluation of  normalized Legendre functions
c                   and their derivatives (with scaling as in ylgndr2s).
c                   Requires prior call to ylgndrfwini.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       (not scaled by sqrt(2*n+1))
c
c       ylgndrufw - faster evaluation of normalized Legendre functions.
c                 Requires prior call to ylgndrfwini.
c
c       ylgndru2fw - faster evaluation of normalized Legendre functions
c                  and their derivatives.
c                  Requires prior call to ylgndrfwini.
c
c       ylgndru2sfw - faster evaluation of  normalized Legendre functions
c                   and their derivatives (with scaling as in ylgndr2s).
c                   Requires prior call to ylgndrfwini.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     CODES
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine ylgndr(nmax, x, y)
      implicit none
c
c     Evaluate normalized Legendre functions defined as:
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                 must be non-negative
c     x                    -1 <= x <= 1
c     y(0:nmax,0:nmax)     resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x)
c     for 0 <= n <= nmax  and  0 <= m <= n.  Other elements of y
c     will contain undefined values.
c
cf2py intent(in) nmax
cf2py intent(in) x
cf2py intent(out) y

      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), u
c
c
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndr2(nmax, x, y, d)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x
cf2py intent(out) y,d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
c
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndr2s(nmax, x, y, d)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     For m>0,
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x
cf2py intent(out) y,d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
      do n=m+2, nmax
        y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
        d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+(1-x**2)*y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+(1-x**2)*y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
      subroutine ylgndru(nmax, x, y)
      implicit none
c
c     Evaluate normalized Legendre functions defined as:
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                 must be non-negative
c     x                    -1 <= x <= 1
c     y(0:nmax,0:nmax)     resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x)
c     for 0 <= n <= nmax  and  0 <= m <= n.  Other elements of y
c     will contain undefined values.
c
cf2py intent(in) nmax
cf2py intent(in) x
cf2py intent(out) y

      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), u
c
c
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c      do n=0, nmax
c	 do m=0, n
c	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
c         enddo
c      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndru2(nmax, x, y, d)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x
cf2py intent(out) y,d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
c
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c      do n=0, nmax
c	 do m=0, n
c	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
c	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
c         enddo
c      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndru2s(nmax, x, y, d)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     For m>0,
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x
cf2py intent(out) y,d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, u2
c
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
      do n=m+2, nmax
        y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
        d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+u2*y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c      do n=0, nmax
c	 do m=0, n
c	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
c	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
c         enddo
c      enddo
      return
      end
c
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
c
      subroutine ylgndr2sm(nmax, x, y, d)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     For m>0,
c          the functions are scaled by 1/sqrt(1-x**2)**m
c          the derivatives are scaled by 1/sqrt(1-x**2)**(m-2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x
cf2py intent(out) y,d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, u2
c
      u=-sqrt((1-x))*sqrt((1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
      do n=m+2, nmax
        y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
        d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+u2*y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       faster version for real argument
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
      subroutine ylgndrini(nmax, rat1, rat2)
      implicit none
c
c     Precompute the recurrence coefficients for the fast
c     evaluation of normalized Legendre functions and their derivatives
c
c     Parameters:
c       nmax                      must be non-negative
c       rat1(0:nmax,0:nmax)       recurrence coefficient
c       rat2(0:nmax,0:nmax)       recurrence coefficient
c
cf2py intent(in) nmax
cf2py intent(out) rat1, rat2
c
      integer nmax, m, n
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      rat1(0,0)=1
      rat2(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  rat1(m,m)=sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  rat2(m,m)=1
	 if (m.lt.nmax)  rat1(m+1,m)=sqrt(2*m+1.0d0)
	 if (m.lt.nmax)  rat2(m+1,m)=1
	 do n=m+2, nmax
	    rat1(n,m)=(2*n-1)
            rat2(n,m)=sqrt((n+m-1.0d0)*(n-m-1.0d0))
	    rat1(n,m)=rat1(n,m)/sqrt(dble(n-m)*(n+m))
	    rat2(n,m)=rat2(n,m)/sqrt(dble(n-m)*(n+m))
         enddo
      enddo
c
c      do m=0, nmax
c	 do n=m, nmax
c            rat1(m,n)=rat1(n,m)
c            rat2(m,n)=rat2(n,m)
c         enddo
c      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndrf(nmax, x, y, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nmax, x, rat1, rat2
cf2py intent(out) y
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2f(nmax, x, y, d, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
c
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2sf(nmax, x, y, d, rat1, rat2)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
      do n=m+2, nmax
        y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndru2sf(nmax, x, y, d, rat1, rat2)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) =  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx =  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, x, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, n, m
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
      do n=m+2, nmax
        y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
c      do n=0, nmax
c	 do m=0, n
c	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
c	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
c        enddo
c      enddo
      return
      end
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       complex valued Legendre functions
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
      subroutine zylgndr(nmax, z, y)
      implicit none
c
c     Evaluate normalized Legendre function for complex argument
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex*16
c     y(0:nmax,0:nmax)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(z)
c     for 0 <= n <= nmax  and  0 <= m <= n.  Other elements of y
c     will contain undefined values.
c
cf2py intent(in) nmax, z
cf2py intent(out) y
c
      integer nmax
      complex*16 z, y(0:nmax,0:nmax)
c
      integer m,n
      complex*16 u
c
      u=-sqrt(1-z*z)
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*z*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine zylgndr2(nmax, z, y, d)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c       for complex argument
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex *16
c     y(0:nmax,0:nmax)      resulting function values, complex *16
c     d(0:nmax,0:nmax)      resulting derivative values, complex *16
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention is valid for the derivative.
c
cf2py intent(in) nmax, z
cf2py intent(out) y, d
c
      integer nmax, m, n
      complex *16 z, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
      u=-sqrt(1-z*z)
      du=-z/u
      y(0,0)=1
      d(0,0)=0
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*z/u**2
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)  d(m+1,m)=(z*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*z*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(z*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine zylgndr2s(nmax, z, y, d)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(z) with m>0
c          the functions are scaled by 1/sqrt(1-z**2)
c          the derivatives are scaled by sqrt(1-z**2)
c
c
c      Ynm(z) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(z)
c
c      d Ynm(z) / dz = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(z) / dz
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex *16
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(z) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, z
cf2py intent(out) y, d
c
      integer nmax, m, n
      complex *16 z, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u
      u=-sqrt(1-z*z)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*sqrt(2*m+1.0d0)
      if (m.lt.nmax)  d(m+1,m)=(z*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
      do n=m+2, nmax
        y(n,m)=((2*n-1)*z*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
        d(n,m)=((2*n-1)*(z*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*z
c
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)
     $      d(m+1,m)=(z*d(m,m)+(1-z**2)*y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*z*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(z*d(n-1,m)+(1-z**2)*y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine zylgndrf(nmax, z, y, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(z) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(z)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex *16
c     y(0:nmax,0:nmax)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(z) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nmax, z, rat1, rat2
cf2py intent(out) y
c
      integer nmax, m, n
      complex *16 z, y(0:nmax,0:nmax), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt(1-z*z)
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*z*y(n-1,m)-rat2(n,m)*y(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine zylgndr2f(nmax, z, y, d, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(z) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(z)
c
c      d Ynm(z) / dz = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(z) / dz
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex *16
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(z) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, z, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, n
      complex *16 z, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt(1-z*z)
      du=z/sqrt(1-z*z)
      y(0,0)=1
      d(0,0)=0
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=(d(m-1,m-1)*u+du)*rat1(m,m)
ccc	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*z/u**2
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)  d(m+1,m)=(z*d(m,m)+y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*z*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(z*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine zylgndr2sf(nmax, z, y, d, rat1, rat2)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(z) with m>0
c          the functions are scaled by 1/sqrt(1-z**2)
c          the derivatives are scaled by sqrt(1-z**2)
c
c
c      Ynm(z) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(z)
c
c      d Ynm(z) / dz = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(z) / dz
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex *16
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(z) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, z, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, n
      complex *16 z, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt(1-z*z)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*rat1(m+1,m)
      if (m.lt.nmax)  d(m+1,m)=(z*d(m,m)+y(m,m))*rat1(m+1,m)
      do n=m+2, nmax
        y(n,m)=rat1(n,m)*z*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(z*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nmax
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*z
c
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)
     $      d(m+1,m)=(z*d(m,m)+(1-z**2)*y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*z*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(z*d(n-1,m)+(1-z**2)*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       complex valued Legendre functions, modified branch cut
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine zylgndrbr(nmax, z, y)
      implicit none
c
c     Evaluate normalized Legendre function for complex argument
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     z                     complex*16
c     y(0:nmax,0:nmax)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(z)
c     for 0 <= n <= nmax  and  0 <= m <= n.  Other elements of y
c     will contain undefined values.
c
c       Modified branch cut at (0,+i)
c
cf2py intent(in) nmax, z
cf2py intent(out) y
c
      integer nmax
      complex*16 z, y(0:nmax,0:nmax)
c
      integer m,n
      complex*16 u
c
      u=-sqrt(1-z*z)
c
c     branch cut at (0,+i), select the lower branch
c     of complex square root
c
      if( imag(1-z*z) .gt. 0 .and. real(1-z*z) .lt. 0) u=+sqrt(1-z*z)
ccc      call prin2('in zylgndrbr, u=*', -u, 2)
ccc      call prin2('in zylgndrbr, 1-z^2=*', 1-z*z, 2)
c
      y(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.lt.nmax)  y(m+1,m)=z*y(m,m)*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*z*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       complex valued Legendre functions (scaled)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
c
      subroutine zylgndrsc(nmax, z,scale, ysc)
      implicit none
c
c     Evaluate scaled versions of zylgndr.  zylgndr is the complex version
c     of the normalized Legendre function.  Scaling removes the possible
c     conditioning errors from zylgndr evaluated at large arguments.
c
c     y_nm = scale**(-n) * ysc_nm
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                      must be non-negative
c     z (complex*16)            complex *16 argument
c     scale (real*8)
c          looks like scale wants to be O(1/z) and scale<1.
c     ysc(0:nmax,0:nmax)        resulting function values
c
cf2py intent(in) nmax, z, scale
cf2py intent(out) ysc
c
      integer nmax
      complex*16 z, ysc(0:nmax,0:nmax)
      real*8 scale
c
      integer m,n
      complex*16 u
c
      u=-sqrt(1-z*z)
      ysc(0,0)=1
      do m=0, nmax
	 if (m.gt.0)  then
            ysc(m,m)=ysc(m-1,m-1)*scale*u*sqrt((2*m-1.0d0)/(2*m))
c           call prinf('m=*',m,1)
c           call prin2('ysc(m,m)=*',ysc(m,m),2)
         endif
	 if (m.lt.nmax)  then
            ysc(m+1,m)=z*scale*ysc(m,m)*sqrt(2*m+1.0d0)
         endif
	 do n=m+2, nmax
	    ysc(n,m)=((2*n-1)*scale*z*ysc(n-1,m) -
     1           sqrt((n+m-1.0d0)*(n-m-1.0d0))*scale**2*ysc(n-2,m))
     2           /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, n
	    ysc(n,m)=ysc(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       truncated recurrences for Legendre functions:
c
c       Ynm(theta) for n = 0,nmax  but m = -m2,...,m2  with m2 < nmax.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c
      subroutine ylgndr2s_trunc(nmax, m2, x, y, d)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, m2, x
cf2py intent(out) y, d
c
      integer nmax, m, m2, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*sqrt(2*m+1.0d0)
      do n=m+2, nmax
        y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
        d(n,m)=((2*n-1)*(x*d(n-1,m)+y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, m2
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*sqrt((2*m-1.0d0)/(2*m))
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*sqrt(2*m+1.0d0)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+(1-x**2)*y(m,m))*sqrt(2*m+1.0d0)
	 do n=m+2, nmax
	    y(n,m)=((2*n-1)*x*y(n-1,m) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*y(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
	    d(n,m)=((2*n-1)*(x*d(n-1,m)+(1-x**2)*y(n-1,m)) -
     1               sqrt((n+m-1.0d0)*(n-m-1.0d0))*d(n-2,m))
     2               /sqrt((n-m+0.0d0)*(n+m))
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, min(n,m2)
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
c
      subroutine ylgndrf_trunc(nmax, m2, x, y, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., min(m2,n).
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nmax, m2, x, rat1, rat2
cf2py intent(out) y
c
      integer nmax, m, m2, n
      real *8 x, y(0:nmax,0:nmax), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do m=0, m2
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, min(n,m2)
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2f_trunc(nmax, m2, x, y, d, rat1, rat2)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., min(m2,n).
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, m2, x, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, m2, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, du
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do m=0, m2
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, min(n,m2)
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2sf_trunc(nmax, m2, x, y, d, rat1, rat2)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nmax
c     and  m = 0, 1,..., min(n,m2).
c
c     Parameters:
c     nmax                  must be non-negative
c     x                     -1 <= x <= 1
c     y(0:nmax,0:nmax)      resulting function values
c     d(0:nmax,0:nmax)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nmax and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nmax, m2, x, rat1, rat2
cf2py intent(out) y, d
c
      integer nmax, m, m2, n
      real *8 x, y(0:nmax,0:nmax), d(0:nmax,0:nmax), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
c
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      if (m.lt.nmax)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
      do n=m+2, nmax
        y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, m2
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nmax)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nmax)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
	 do n=m+2, nmax
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nmax
	 do m=0, min(n,m2)
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       Symmetries for Legendre functions
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
        subroutine ylgndrpm(nterms,y)
        implicit none
        integer n,m,nterms
        real *8 y(0:nterms,0:nterms)
c
c       Given Y_nm(x), return Y_nm(-x)
c
        do n=0,nterms
           do m=0,n
              if( mod(n+m,2) .eq. 1 ) y(n,m)=-y(n,m)
           enddo
        enddo
c
        return
        end
c
c
c
c
c
        subroutine ylgndr2pm(nterms,y,d)
        implicit none
        integer nterms,n,m
        real *8 y(0:nterms,0:nterms)
        real *8 d(0:nterms,0:nterms)
c
c       Given Y_nm(x), return Y_nm(-x)
c       Given Y'_nm(x), return Y'_nm(-x)
c
cf2py intent(in) nterms, y
cf2py intent(out) d
c
        do n=0,nterms
           do m=0,n
              if( mod(n+m,2) .eq. 1 ) y(n,m)=-y(n,m)
              if( mod(n+m,2) .eq. 0 ) d(n,m)=-d(n,m)
           enddo
        enddo
c
        return
        end
c
c
c
c
c
        subroutine ylgndrpm_opt(nterms,y)
        implicit none
        integer nterms,n,m
        real *8 y(0:nterms,0:nterms)
c
c       Given Y_nm(x), return Y_nm(-x)
c
        do n=0,nterms,2
           do m=1,n,2
              y(n,m)=-y(n,m)
           enddo
        enddo
c
        do n=1,nterms,2
           do m=0,n,2
              y(n,m)=-y(n,m)
           enddo
        enddo
c
        return
        end
c
c
c
c
c
        subroutine ylgndr2pm_opt(nterms,y,d)
        implicit none
        integer nterms,n,m
        real *8 y(0:nterms,0:nterms)
        real *8 d(0:nterms,0:nterms)
c
c       Given Y_nm(x), return Y_nm(-x)
c       Given Y'_nm(x), return Y'_nm(-x)
c
cf2py intent(in) nterms, y
cf2py intent(out) d
c
        do n=0,nterms,2
           do m=0,n,2
              d(n,m)=-d(n,m)
           enddo
           do m=1,n,2
              y(n,m)=-y(n,m)
           enddo
        enddo
c
        do n=1,nterms,2
           do m=0,n,2
              y(n,m)=-y(n,m)
           enddo
           do m=1,n,2
              d(n,m)=-d(n,m)
           enddo
        enddo
c
        return
        end
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       fast version of real valued Legendre functions, with storage
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
      subroutine ylgndrfwini(nmax, w, lw, lused)
      implicit none
c
c     Precompute the recurrence coefficients for the fast
c     evaluation of normalized Legendre functions and their derivatives
c
c     Parameters:
c       nmax             must be non-negative
c       w                  contains rat1 and rat2 arrays
c
      integer nmax,irat1,irat2,lw,lused
      real *8 w(*)
c
cf2py intent(in) nmax, lw
cf2py intent(out) w, lused
c
      irat1=1
      irat2=1+(nmax+1)**2
      lused=2*(nmax+1)**2
      if( lused .gt. lw ) return

      call ylgndrini(nmax, w(irat1), w(irat2))
      return
      end
c
c
c
c
c
      subroutine ylgndrfw(nterms, x, y, w, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)          resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndrfw0(nterms, x, y, w(irat1), w(irat2), nmax)
      else
        call ylgndr(nterms, x, y)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2fw(nterms, x, y, d, w, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y, d
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndr2fw0(nterms, x, y, d, w(irat1), w(irat2), nmax)
      else
        call ylgndr2(nterms, x, y, d)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2sfw(nterms, x, y, d, w, nmax)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y, d
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndr2sfw0(nterms, x, y, d, w(irat1), w(irat2), nmax)
      else
        call ylgndr2s(nterms, x, y, d)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndrfw0(nterms, x, y, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y
c
      integer nterms,nmax,m,n
      real *8 x, y(0:nterms,0:nterms), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do m=0, nterms
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 do n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         enddo
      enddo
c
      do n=0, nterms
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2fw0(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, m, n
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, du
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do m=0, nterms
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nterms)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
	 do n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      do n=0, nterms
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndr2sfw0(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt(2n+1)  sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, m, n
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      if (m.lt.nterms)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
      do n=m+2, nterms
        y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
      enddo
c
c       ... then, evaluate scaled associated Legendre functions
c
      do m=1, nterms
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nterms)
ccc     $      d(m+1,m)=(x*d(m,m)+(1-x**2)*y(m,m))*rat1(m+1,m)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
	 do n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
ccc	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+(1-x**2)*y(n-1,m))-
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
         enddo
      enddo
      do n=0, nterms
	 do m=0, n
	    y(n,m)=y(n,m)*sqrt(2*n+1.0d0)
	    d(n,m)=d(n,m)*sqrt(2*n+1.0d0)
         enddo
      enddo
      return
      end
c
c
c
c
c
      subroutine ylgndrufw(nterms, x, y, w, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)          resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndrufw0(nterms, x, y, w(irat1), w(irat2), nmax)
      else
        call ylgndru(nterms, x, y)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndru2fw(nterms, x, y, d, w, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y, d
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndru2fw0(nterms, x, y, d, w(irat1), w(irat2), nmax)
      else
        call ylgndru2(nterms, x, y, d)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndru2sfw(nterms, x, y, d, w, nmax)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, w, nmax
cf2py intent(out) y, d
c
      integer nterms,nmax,irat1,irat2
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), w(*)
c
      if( nterms .le. nmax ) then
        irat1=1
        irat2=1+(nmax+1)**2
        call ylgndru2sfw0(nterms, x, y, d, w(irat1), w(irat2), nmax)
      else
        call ylgndru2s(nterms, x, y, d)
      endif
c
      return
      end
c
c
c
c
c
      subroutine ylgndrufw0_old(nterms, x, y, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      y(0,0)=1
      do 10 m=0, nterms
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 do 20 n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
 20      continue
 10   continue
c
      return
      end
c
c
c
c
c
      subroutine ylgndrufw0(nterms, x, y, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), u
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)

      y(0,0)=1
      if( nterms .eq. 0 ) return

      y(1,0)=x*y(0,0)*rat1(1,0)

      u=-sqrt((1-x)*(1+x))
      do m=1, nterms-1
	 y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      enddo

      y(nterms,nterms)=y(nterms-1,nterms-1)*u*rat1(nterms,nterms)

c      do m=0, nterms-1
c	 do n=m+2, nterms
c	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
c         enddo
c      enddo

      do n=2, nterms
         do m=0, n-2
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndru2fw0_old(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, du
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      du=x/sqrt((1-x)*(1+x))
      y(0,0)=1
      d(0,0)=0
      do 10 m=0, nterms
	 if (m.gt.0)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x/u**2
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nterms)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
	 do 20 n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
 20      continue
 10   continue
c
      return
      end
c
c
c
c
c
      subroutine ylgndru2fw0(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate normalized Legendre functions and their derivatives
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)

      y(0,0)=1
      d(0,0)=0
      if( nterms .eq. 0 ) return

      y(1,0)=x*y(0,0)*rat1(1,0)
      d(1,0)=(x*d(0,0)+y(0,0))*rat1(1,0)

      u=-sqrt((1-x)*(1+x))
      u2=(1-x)*(1+x)
      do m=1, nterms-1
         y(m,m)=y(m-1,m-1)*u*rat1(m,m)
         d(m,m)=y(m,m)*(-m)*x/u2
         y(m+1,m)=x*y(m,m)*rat1(m+1,m)
         d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
      enddo

      y(nterms,nterms)=y(nterms-1,nterms-1)*u*rat1(nterms,nterms)
      d(nterms,nterms)=y(nterms,nterms)*(-nterms*x)

      do n=2, nterms
         do m=0, n-2
         y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      return
      end
c
c
c
c
c
      subroutine ylgndru2sfw0_old(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)
      u=-sqrt((1-x)*(1+x))
      u2 = (1-x)*(1+x)
      y(0,0)=1
      d(0,0)=0
c
c       ... first, evaluate standard Legendre polynomials
c
      m=0
      if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
      if (m.lt.nterms)  d(m+1,m)=(x*d(m,m)+y(m,m))*rat1(m+1,m)
      do 120 n=m+2, nterms
        y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
        d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
120   continue
c
c       ... then, evaluate scaled associated Legendre functions
c
      do 210 m=1, nterms
c
	 if (m.eq.1)  y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
	 if (m.gt.1)  y(m,m)=y(m-1,m-1)*u*rat1(m,m)
	 if (m.gt.0)  d(m,m)=y(m,m)*(-m)*x
c
	 if (m.lt.nterms)  y(m+1,m)=x*y(m,m)*rat1(m+1,m)
	 if (m.lt.nterms)
     $      d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
	 do 220 n=m+2, nterms
	    y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
	    d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-
     $         rat2(n,m)*d(n-2,m)
220      continue
210   continue
c
      return
      end
c
c
c

      subroutine ylgndru2sfw0(nterms, x, y, d, rat1, rat2, nmax)
      implicit none
c
c     Evaluate scaled normalized Legendre functions and their derivatives
c
c     Only for Ynm(x) with m>0
c          the functions are scaled by 1/sqrt(1-x**2)
c          the derivatives are scaled by sqrt(1-x**2)
c
c
c      Ynm(x) = sqrt( (n-m)!/ (n+m)! ) Pnm(x)
c
c      d Ynm(x) / dx = sqrt( (n-m)!/ (n+m)! ) d Pnm(x) / dx
c
c     for n = 0, 1, 2,..., nterms
c     and  m = 0, 1,..., n.
c
c     Parameters:
c     nterms                      must be non-negative
c     x                         -1 <= x <= 1
c     y(0:nterms,0:nterms)      resulting function values
c     d(0:nterms,0:nterms)      resulting derivative values
c
c     Upon return, y(n,m) will contain the function value Ynm(x) for 0
c     <= n <= nterms and 0 <= m <= n.  Other elements of y will contain
c     undefined values. The same convention for the derivatives.
c
cf2py intent(in) nterms, x, rat1, rat2, nmax
cf2py intent(out) y, d
c
      integer nterms, nmax, n, m
      real *8 x, y(0:nterms,0:nterms), d(0:nterms,0:nterms), u, u2
      real *8 rat1(0:nmax,0:nmax), rat2(0:nmax,0:nmax)

      y(0,0)=1
      d(0,0)=0
      if( nterms .eq. 0 ) return

      y(1,0)=x*y(0,0)*rat1(1,0)
      d(1,0)=(x*d(0,0)+y(0,0))*rat1(1,0)

      u=-sqrt((1-x)*(1+x))
      u2=(1-x)*(1+x)
      do m=1, nterms-1
         if( m .eq. 1 ) y(m,m)=y(m-1,m-1)*(-1)*rat1(m,m)
         if( m .gt. 1 ) y(m,m)=y(m-1,m-1)*u*rat1(m,m)
         d(m,m)=y(m,m)*(-m)*x
         y(m+1,m)=x*y(m,m)*rat1(m+1,m)
         d(m+1,m)=(x*d(m,m)+u2*y(m,m))*rat1(m+1,m)
      enddo

      y(nterms,nterms)=y(nterms-1,nterms-1)*u*rat1(nterms,nterms)
      d(nterms,nterms)=y(nterms,nterms)*(-nterms*x)

      do n=2, nterms
         m=0
         y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         d(n,m)=rat1(n,m)*(x*d(n-1,m)+y(n-1,m))-rat2(n,m)*d(n-2,m)
         do m=1, n-2
         y(n,m)=rat1(n,m)*x*y(n-1,m)-rat2(n,m)*y(n-2,m)
         d(n,m)=rat1(n,m)*(x*d(n-1,m)+u2*y(n-1,m))-rat2(n,m)*d(n-2,m)
         enddo
      enddo
c
      return
      end