rmn2l Subroutine

public subroutine rmn2l(m, n, c, x, Df, Kd, R2f, R2d, Id)

Compute prolate and oblate spheroidal radial functions of the second kind for given m, n, c and a large cx

Arguments

Type IntentOptional Attributes Name
integer :: m
integer :: n
real(kind=wp) :: c
real(kind=wp) :: x
real(kind=wp), dimension(200) :: Df
integer :: Kd
real(kind=wp) :: R2f
real(kind=wp) :: R2d
integer :: Id

Calls

proc~~rmn2l~~CallsGraph proc~rmn2l specfun_module::rmn2l proc~sphy specfun_module::sphy proc~rmn2l->proc~sphy

Called by

proc~~rmn2l~~CalledByGraph proc~rmn2l specfun_module::rmn2l proc~rswfo specfun_module::rswfo proc~rswfo->proc~rmn2l proc~rswfp specfun_module::rswfp proc~rswfp->proc~rmn2l

Source Code

   subroutine rmn2l(m,n,c,x,Df,Kd,R2f,R2d,Id)

      real(wp) :: a0 , b0 , c , cx , Df , dy , eps1 , eps2 , &
                  r , r0 , R2d , R2f , reg , suc , sud , sw , sy , &
                  x
      integer :: Id , id1 , id2 , ip , j , k , Kd , l , lg , m , n , nm , &
                 nm1 , nm2 , np
      dimension Df(200) , sy(0:251) , dy(0:251)

      real(wp),parameter :: eps = 1.0e-14_wp

      ip = 1
      nm1 = int((n-m)/2)
      if ( n-m==2*nm1 ) ip = 0
      nm = 25 + nm1 + int(c)
      reg = 1.0_wp
      if ( m+nm>80 ) reg = 1.0e-200_wp
      nm2 = 2*nm + m
      cx = c*x
      call sphy(nm2,cx,nm2,sy,dy)
      r0 = reg
      do j = 1 , 2*m + ip
         r0 = r0*j
      enddo
      r = r0
      suc = r*Df(1)
      sw = 0.0_wp
      do k = 2 , nm
         r = r*(m+k-1.0_wp)*(m+k+ip-1.5_wp)/(k-1.0_wp)/(k+ip-1.5_wp)
         suc = suc + r*Df(k)
         if ( k>nm1 .and. abs(suc-sw)<abs(suc)*eps ) exit
         sw = suc
      enddo
      a0 = (1.0_wp-Kd/(x*x))**(0.5_wp*m)/suc
      R2f = 0.0_wp
      eps1 = 0.0_wp
      np = 0
      do k = 1 , nm
         l = 2*k + m - n - 2 + ip
         lg = 1
         if ( l/=4*int(l/4) ) lg = -1
         if ( k==1 ) then
            r = r0
         else
            r = r*(m+k-1.0_wp)*(m+k+ip-1.5_wp)/(k-1.0_wp)/(k+ip-1.5_wp)
         endif
         np = m + 2*k - 2 + ip
         R2f = R2f + lg*r*(Df(k)*sy(np))
         eps1 = abs(R2f-sw)
         if ( k>nm1 .and. eps1<abs(R2f)*eps ) exit
         sw = R2f
      enddo
      id1 = int(log10(eps1/abs(R2f)+eps))
      R2f = R2f*a0
      if ( np>=nm2 ) then
         Id = 10
         return
      endif
      b0 = Kd*m/x**3.0_wp/(1.0-Kd/(x*x))*R2f
      sud = 0.0_wp
      eps2 = 0.0_wp
      do k = 1 , nm
         l = 2*k + m - n - 2 + ip
         lg = 1
         if ( l/=4*int(l/4) ) lg = -1
         if ( k==1 ) then
            r = r0
         else
            r = r*(m+k-1.0_wp)*(m+k+ip-1.5_wp)/(k-1.0_wp)/(k+ip-1.5_wp)
         endif
         np = m + 2*k - 2 + ip
         sud = sud + lg*r*(Df(k)*dy(np))
         eps2 = abs(sud-sw)
         if ( k>nm1 .and. eps2<abs(sud)*eps ) exit
         sw = sud
      enddo
      R2d = b0 + a0*c*sud
      id2 = int(log10(eps2/abs(sud)+eps))
      Id = max(id1,id2)

   end subroutine rmn2l