cva1 Subroutine

public subroutine cva1(Kd, m, q, Cv)

Compute a sequence of characteristic values of Mathieu functions

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Kd

Case code:

  • KD=1 for cem(x,q) ( m = 0,2,4,… )
  • KD=2 for cem(x,q) ( m = 1,3,5,… )
  • KD=3 for sem(x,q) ( m = 1,3,5,… )
  • KD=4 for sem(x,q) ( m = 2,4,6,… )
integer, intent(in) :: m

Maximum order of Mathieu functions
real(kind=wp), intent(in) :: q

Parameter of Mathieu functions
real(kind=wp), intent(out) :: Cv(200)

CV(I) --- Characteristic values; I = 1,2,3,...

  • For KD=1, CV(1), CV(2), CV(3),..., correspond to the characteristic values of cem for m = 0,2,4,...
  • For KD=2, CV(1), CV(2), CV(3),..., correspond to the characteristic values of cem for m = 1,3,5,...
  • For KD=3, CV(1), CV(2), CV(3),..., correspond to the characteristic values of sem for m = 1,3,5,...
  • For KD=4, CV(1), CV(2), CV(3),..., correspond to the characteristic values of sem for m = 0,2,4,...

Source Code

   subroutine cva1(Kd,m,q,Cv)

      integer,intent(in) :: Kd !! Case code:
                               !!
                               !!  * `KD=1` for `cem(x,q)`  ( `m = 0,2,4,…` )
                               !!  * `KD=2` for `cem(x,q)`  ( `m = 1,3,5,…` )
                               !!  * `KD=3` for `sem(x,q)`  ( `m = 1,3,5,…` )
                               !!  * `KD=4` for `sem(x,q)`  ( `m = 2,4,6,…` )
      integer,intent(in) :: m !! Maximum order of Mathieu functions
      real(wp),intent(in) :: q !! Parameter of Mathieu functions
      real(wp),intent(out) :: Cv(200) !! CV(I) --- Characteristic values; I = 1,2,3,...
                                      !!
                                      !! * For `KD=1, CV(1), CV(2), CV(3),...`, correspond to
                                      !!   the characteristic values of `cem` for `m = 0,2,4,...`
                                      !! * For `KD=2, CV(1), CV(2), CV(3),...`, correspond to
                                      !!   the characteristic values of `cem` for `m = 1,3,5,...`
                                      !! * For `KD=3, CV(1), CV(2), CV(3),...`, correspond to
                                      !!   the characteristic values of `sem` for `m = 1,3,5,...`
                                      !! * For `KD=4, CV(1), CV(2), CV(3),...`, correspond to
                                      !!   the characteristic values of `sem` for `m = 0,2,4,...`

      real(wp),dimension(500) :: d , e , f
      real(wp),dimension(200) :: g , h
      real(wp) :: s , t , t1 , x1 , xa , xb
      integer :: i , ic , icm , j , k , k1 , nm , nm1

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

      icm = int(m/2) + 1
      if ( Kd==4 ) icm = m/2
      if ( q/=0.0_wp ) then
         nm = int(10.0_wp+1.5_wp*m+0.5_wp*q)
         e(1) = 0.0_wp
         f(1) = 0.0_wp
         if ( Kd==1 ) then
            d(1) = 0.0_wp
            do i = 2 , nm
               d(i) = 4.0_wp*(i-1.0_wp)**2
               e(i) = q
               f(i) = q*q
            enddo
            e(2) = sq2*q
            f(2) = 2.0_wp*q*q
         elseif ( Kd/=4 ) then
            d(1) = 1.0_wp + (-1)**Kd*q
            do i = 2 , nm
               d(i) = (2.0_wp*i-1.0_wp)**2
               e(i) = q
               f(i) = q*q
            enddo
         else
            d(1) = 4.0_wp
            do i = 2 , nm
               d(i) = 4.0_wp*i*i
               e(i) = q
               f(i) = q*q
            enddo
         endif
         xa = d(nm) + abs(e(nm))
         xb = d(nm) - abs(e(nm))
         nm1 = nm - 1
         do i = 1 , nm1
            t = abs(e(i)) + abs(e(i+1))
            t1 = d(i) + t
            if ( xa<t1 ) xa = t1
            t1 = d(i) - t
            if ( t1<xb ) xb = t1
         enddo
         do i = 1 , icm
            g(i) = xa
            h(i) = xb
         enddo
         do k = 1 , icm
            do k1 = k , icm
               if ( g(k1)<g(k) ) then
                  g(k) = g(k1)
                  exit
               endif
            enddo
            if ( k/=1 .and. h(k)<h(k-1) ) h(k) = h(k-1)
            do
               x1 = (g(k)+h(k))/2.0_wp
               Cv(k) = x1
               if ( abs((g(k)-h(k))/x1)<eps ) then
                  Cv(k) = x1
                  exit
               else
                  j = 0
                  s = 1.0_wp
                  do i = 1 , nm
                     if ( s==0.0_wp ) s = s + 1.0e-30_wp
                     t = f(i)/s
                     s = d(i) - t - x1
                     if ( s<0.0_wp ) j = j + 1
                  enddo
                  if ( j<k ) then
                     h(k) = x1
                  else
                     g(k) = x1
                     if ( j>=icm ) then
                        g(icm) = x1
                     else
                        if ( h(j+1)<x1 ) h(j+1) = x1
                        if ( x1<g(j) ) g(j) = x1
                     endif
                  endif
               endif
            end do
         enddo
      elseif ( Kd==1 ) then
         do ic = 1 , icm
            Cv(ic) = 4.0_wp*(ic-1.0_wp)**2
         enddo
      elseif ( Kd/=4 ) then
         do ic = 1 , icm
            Cv(ic) = (2.0_wp*ic-1.0_wp)**2
         enddo
      else
         do ic = 1 , icm
            Cv(ic) = 4.0_wp*ic*ic
         enddo
      endif

   end subroutine cva1