cbk Subroutine

public subroutine cbk(m, n, c, Cv, Qt, Ck, Bk)

Compute coefficient Bk's for oblate radial functions with a small argument

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: m
integer, intent(in) :: n
real(kind=wp), intent(in) :: c
real(kind=wp), intent(in) :: Cv
real(kind=wp), intent(in) :: Qt
real(kind=wp), intent(in) :: Ck(200)
real(kind=wp), intent(out) :: Bk(200)

Called by

proc~~cbk~~CalledByGraph proc~cbk specfun_module::cbk proc~rmn2so specfun_module::rmn2so proc~rmn2so->proc~cbk proc~rswfo specfun_module::rswfo proc~rswfo->proc~rmn2so

Source Code

   subroutine cbk(m,n,c,Cv,Qt,Ck,Bk)

      integer,intent(in) :: m
      integer,intent(in) :: n
      real(wp),intent(in) :: c
      real(wp),intent(in) :: Cv
      real(wp),intent(in) :: Qt
      real(wp),intent(in) :: Ck(200)
      real(wp),intent(out) :: Bk(200)

      real(wp) :: r1 , s1 , sw , t , u(200) , v(200) , w(200)
      integer :: i , i1 , ip , j , k , n2 , nm

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

      ip = 1
      if ( n-m==2*int((n-m)/2) ) ip = 0
      nm = 25 + int(0.5_wp*(n-m)+c)
      u(1) = 0.0_wp
      n2 = nm - 2
      do j = 2 , n2
         u(j) = c*c
      enddo
      do j = 1 , n2
         v(j) = (2.0_wp*j-1.0_wp-ip)*(2.0_wp*(j-m)-ip) + m*(m-1.0_wp) - Cv
      enddo
      do j = 1 , nm - 1
         w(j) = (2.0_wp*j-ip)*(2.0_wp*j+1.0_wp-ip)
      enddo
      if ( ip==0 ) then
         sw = 0.0_wp
         do k = 0 , n2 - 1
            s1 = 0.0_wp
            i1 = k - m + 1
            do i = i1 , nm
               if ( i>=0 ) then
                  r1 = 1.0_wp
                  do j = 1 , k
                     r1 = r1*(i+m-j)/j
                  enddo
                  s1 = s1 + Ck(i+1)*(2.0_wp*i+m)*r1
                  if ( abs(s1-sw)<abs(s1)*eps ) exit
                  sw = s1
               endif
            enddo
            Bk(k+1) = Qt*s1
         enddo
      elseif ( ip==1 ) then
         sw = 0.0_wp
         do k = 0 , n2 - 1
            s1 = 0.0_wp
            i1 = k - m + 1
            do i = i1 , nm
               if ( i>=0 ) then
                  r1 = 1.0_wp
                  do j = 1 , k
                     r1 = r1*(i+m-j)/j
                  enddo
                  if ( i>0 ) s1 = s1 + Ck(i)*(2.0_wp*i+m-1)*r1
                  s1 = s1 - Ck(i+1)*(2.0_wp*i+m)*r1
                  if ( abs(s1-sw)<abs(s1)*eps ) exit
                  sw = s1
               endif
            enddo
            Bk(k+1) = Qt*s1
         enddo
      endif
      w(1) = w(1)/v(1)
      Bk(1) = Bk(1)/v(1)
      do k = 2 , n2
         t = v(k) - w(k-1)*u(k)
         w(k) = w(k)/t
         Bk(k) = (Bk(k)-Bk(k-1)*u(k))/t
      enddo
      do k = n2 - 1 , 1 , -1
         Bk(k) = Bk(k) - w(k)*Bk(k+1)
      enddo

   end subroutine cbk