sdmn Subroutine

subroutine sdmn ( m, n, c, cv, kd, df )

  !*****************************************************************************80
  !
  !! SDMN: expansion coefficients for prolate and oblate spheroidal functions.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    29 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) M, the mode parameter.
  !
  !    Input, integer ( kind = 4 ) N, the mode parameter.
  !
  !    Input, real ( kind = 8 ) C, the spheroidal parameter.
  !
  !    Input, real ( kind = 8 ) CV, the characteristic value.
  !
  !    Input, integer ( kind = 4 ) KD, the function code.
  !    1, the prolate function.
  !    -1, the oblate function.
  !
  !    Output, real ( kind = 8 ) DF(*), expansion coefficients;
  !    DF(1), DF(2), ... correspond to d0, d2, ... for even n-m and d1,
  !    d3, ... for odd n-m
  !
  implicit none

  real ( kind = 8 ) a(200)
  real ( kind = 8 ) c
  real ( kind = 8 ) cs
  real ( kind = 8 ) cv
  real ( kind = 8 ) d(200)
  real ( kind = 8 ) d2k
  real ( kind = 8 ) df(200)
  real ( kind = 8 ) dk0
  real ( kind = 8 ) dk1
  real ( kind = 8 ) dk2
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  real ( kind = 8 ) fl
  real ( kind = 8 ) fs
  real ( kind = 8 ) g(200)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ip
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k1
  integer ( kind = 4 ) kb
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm
  real ( kind = 8 ) r1
  real ( kind = 8 ) r3
  real ( kind = 8 ) r4
  real ( kind = 8 ) s0
  real ( kind = 8 ) su1
  real ( kind = 8 ) su2
  real ( kind = 8 ) sw

  nm = 25 + int ( 0.5D+00 * ( n - m ) + c )

  if ( c < 1.0D-10 ) then
     do i = 1, nm
        df(i) = 0D+00
     end do
     df((n-m)/2+1) = 1.0D+00
     return
  end if

  cs = c * c * kd

  if ( n - m == 2 * int ( ( n - m ) / 2 ) ) then
     ip = 0
  else
     ip = 1
  end if

  do i = 1, nm + 2
     if ( ip == 0 ) then
        k = 2 * ( i - 1 )
     else
        k = 2 * i - 1
     end if
     dk0 = m + k
     dk1 = m + k + 1
     dk2 = 2 * ( m + k )
     d2k = 2 * m + k
     a(i) = ( d2k + 2.0D+00 ) * ( d2k + 1.0D+00 ) &
          / ( ( dk2 + 3.0D+00 ) * ( dk2 + 5.0D+00 ) ) * cs
     d(i) = dk0 * dk1 &
          + ( 2.0D+00 * dk0 * dk1 - 2.0D+00 * m * m - 1.0D+00 ) &
          / ( ( dk2 - 1.0D+00 ) * ( dk2 + 3.0D+00 ) ) * cs
     g(i) = k * ( k - 1.0D+00 ) / ( ( dk2 - 3.0D+00 ) &
          * ( dk2 - 1.0D+00 ) ) * cs
  end do

  fs = 1.0D+00
  f1 = 0.0D+00
  f0 = 1.0D-100
  kb = 0
  df(nm+1) = 0.0D+00

  do k = nm, 1, -1

     f = - ( ( d(k+1) - cv ) * f0 + a(k+1) * f1 ) / g(k+1)

     if ( abs ( df(k+1) ) < abs ( f ) ) then

        df(k) = f
        f1 = f0
        f0 = f
        if ( 1.0D+100 < abs ( f ) ) then
           do k1 = k, nm
              df(k1) = df(k1) * 1.0D-100
           end do
           f1 = f1 * 1.0D-100
           f0 = f0 * 1.0D-100
        end if

     else

        kb = k
        fl = df(k+1)
        f1 = 1.0D-100
        f2 = - ( d(1) - cv ) / a(1) * f1
        df(1) = f1

        if ( kb == 1 ) then

           fs = f2

        else if ( kb == 2 ) then

           df(2) = f2
           fs = - ( ( d(2) - cv ) * f2 + g(2) * f1 ) / a(2)

        else 

           df(2) = f2
           do j = 3, kb + 1
              f = - ( ( d(j-1) - cv ) * f2 + g(j-1) * f1 ) / a(j-1)
              if ( j <= kb ) then
                 df(j) = f
              end if
              if ( 1.0D+100 < abs ( f ) ) then
                 do k1 = 1, j
                    df(k1) = df(k1) * 1.0D-100
                 end do
                 f = f * 1.0D-100
                 f2 = f2 * 1.0D-100
              end if
              f1 = f2
              f2 = f
           end do
           fs = f

        end if

        exit

     end if

  end do

  su1 = 0.0D+00

  r1 = 1.0D+00
  do j = m + ip + 1, 2 * ( m + ip )
     r1 = r1 * j
  end do

  su1 = df(1) * r1
  do k = 2, kb
     r1 = - r1 * ( k + m + ip - 1.5D+00 ) / ( k - 1.0D+00 )
     su1 = su1 + r1 * df(k)
  end do

  su2 = 0.0D+00
  do k = kb + 1, nm
     if ( k /= 1 ) then
        r1 = - r1 * ( k + m + ip - 1.5D+00 ) / ( k - 1.0D+00 )
     end if
     su2 = su2 + r1 * df(k)
     if ( abs ( sw - su2 ) < abs ( su2 ) * 1.0D-14 ) then
        exit
     end if
     sw = su2
  end do

  r3 = 1.0D+00
  do j = 1, ( m + n + ip ) / 2
     r3 = r3 * ( j + 0.5D+00 * ( n + m + ip ) )
  end do

  r4 = 1.0D+00
  do j = 1, ( n - m - ip ) / 2
     r4 = -4.0D+00 * r4 * j
  end do

  s0 = r3 / ( fl * ( su1 / fs ) + su2 ) / r4
  do k = 1, kb
     df(k) = fl / fs * s0 * df(k)
  end do

  do k = kb + 1, nm
     df(k) = s0 * df(k)
  end do

  return
end subroutine sdmn