program main

c*********************************************************************72
c
cc gegenbauer_polynomial_test() tests gegenbauer_polynomial().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    24 April 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'gegenbauer_polynomial_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test gegenbauer_polynomial().'

      call gegenbauer_alpha_check_test ( )
      call gegenbauer_ek_compute_test ( )
      call gegenbauer_integral_test ( )
      call gegenbauer_polynomial_value_test ( )
      call gegenbauer_ss_compute_test ( )
      call gegenbauer_to_monomial_matrix_test ( )

      call imtqlx_test ( )

      call monomial_to_gegenbauer_matrix_test ( )

      call r8_hyper_2f1_test ( )
      call r8_psi_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'gegenbauer_polynomial_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop 0
      end
      subroutine gegenbauer_alpha_check_test ( )

c*********************************************************************72
c
cc gegenbauer_alpha_check_test() compares gegenbauer_alpha_check().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision alpha
      logical check
      integer n
      double precision r

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_alpha_check_test():'
      write ( *, '(a)' ) 
     &  '  gegenbauer_alpha_check() checks that ALPHA is legal.'
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '       ALPHA   Check?'
      write ( *, '(a)' ) ''

      do n = 1, 10

        call random_number ( harvest = r )
        alpha = -5.0D+00 + r * 10.0D+00
        call gegenbauer_alpha_check ( alpha, check )
        write ( *, '(2x,f10.4,7x,l1)' ) alpha, check

      end do

      return
      end
      subroutine gegenbauer_ek_compute_test ( )

c*********************************************************************72
c
cc gegenbauer_ek_compute_test() tests gegenbauer_ek_compute().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    19 November 2015
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision alpha
      integer i
      integer n
      double precision, allocatable, dimension ( : ) :: w
      double precision, allocatable, dimension ( : ) :: x

      alpha = 0.50D+00

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_ek_compute_test():'
      write ( *, '(a)' ) 
     &  '  gegenbauer_ek_compute() computes a Gauss-Gegenbauer rule;'
      write ( *, '(a)' ) ''
      write ( *, '(a,g14.6)' ) '  Using parameter ALPHA = ', alpha
      write ( *, '(a,g14.6,a,g14.6)' ) 
     &  '  Integration interval is [-1,+1].'
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '                 W                         X'
      write ( *, '(a)' ) ''

      do n = 1, 10

        allocate ( w(n) )
        allocate ( x(n) )

        call gegenbauer_ek_compute ( n, alpha, x, w )
     
        write ( *, '(a)' ) ''
        do i = 1, n
          write ( *, '(10x,2x,g24.16,2x,g24.16)' ) w(i), x(i)
        end do

        deallocate ( w )
        deallocate ( x )

      end do

      return
      end
      subroutine gegenbauer_integral_test ( )

c*********************************************************************72
c
cc gegenbauer_integral_test() tests gegenbauer_integral().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    14 June 2015
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision alpha
      integer n
      double precision value

      alpha = 0.50D+00

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_integral_test():'
      write ( *, '(a)' ) '  gegenbauer_integral() evaluates'
      write ( *, '(a)' ) 
     &  '  Integral ( -1 < x < +1 ) x^n * (1-x^2)^(alpha-1/2) dx'
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '         N         Value'
      write ( *, '(a)' ) ''

      do n = 0, 10
     
        call gegenbauer_integral ( n, alpha, value )

        write ( *, '(2x,i8,2x,g24.16)' ) n, value

      end do
     
      return
      end
      subroutine gegenbauer_polynomial_value_test ( )

c*********************************************************************72
c
cc gegenbauer_polynomial_value_test() tests gegenbauer_polynomial_value().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    23 November 2015
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision alpha
      double precision, allocatable :: c(:,:)
      double precision fx(1)
      integer m
      integer n
      integer n_data
      double precision x(1)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_polynomial_value_test():'
      write ( *, '(a)' ) '  gegenbauer_polynomial_value() evaluates '
      write ( *, '(a)' ) '  a Gegenbauer polynomial.'
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 
     &  '       M     ALPHA         X           GPV    GEGENBAUER'
      write ( *, '(a)' ) ''

      n = 1

      n_data = 0

      do

        call gegenbauer_polynomial_values ( n_data, m, alpha, x, fx )

        if ( n_data == 0 ) then
          exit
        end if

        allocate ( c(0:m,1:n) )

        call gegenbauer_polynomial_value ( m, n, alpha, x, c )

        write ( *, '(2x,i6,2x,f8.2,2x,f8.2,2x,f12.4,2x,f12.4)' ) 
     &    m, alpha, x(1), fx, c(m,1)

        deallocate ( c )

      end do

      return
      end
      subroutine gegenbauer_ss_compute_test ( )

c*********************************************************************72
c
cc gegenbauer_ss_compute_test() tests gegenbauer_ss_compute().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    18 June 2015
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision alpha
      integer i
      integer n
      double precision, allocatable, dimension ( : ) :: w
      double precision, allocatable, dimension ( : ) :: x

      alpha = 0.50D+00

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_ss_compute_test():'
      write ( *, '(a)' ) 
     &  '  gegenbauer_ss_compute() computes a Gauss-Gegenbauer rule;'
      write ( *, '(a)' ) ''
      write ( *, '(a,g14.6)' ) '  Using parameter ALPHA = ', alpha
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '                 W                         X'
      write ( *, '(a)' ) ''

      do n = 1, 10

        allocate ( w(n) )
        allocate ( x(n) )

        call gegenbauer_ss_compute ( n, alpha, x, w )
     
        write ( *, '(a)' ) ''
        do i = 1, n
          write ( *, '(10x,2x,g24.16,2x,g24.16)' ) w(i), x(i)
        end do

        deallocate ( w )
        deallocate ( x )

      end do

      return
      end
      subroutine gegenbauer_to_monomial_matrix_test ( )

c*********************************************************************72
c
cc gegenbauer_to_monomial_matrix_test() tests gegenbauer_to_monomial_matrix().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 5

      double precision A(n,n)
      double precision alpha
      double precision g(n)
      integer i
      double precision mono(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'gegenbauer_to_monomial_matrix_test():'
      write ( *, '(a)' ) 
     &  '  gegenbauer_to_monomial_matrix() evaluates the matrix'
      write ( *, '(a)' ) 
     &  '  which converts Gegenbauer polyjomial coefficients'
      write ( *, '(a)' ) '  to monomial coefficients.'
      
      alpha = 0.5D+00

      call gegenbauer_to_monomial_matrix ( n, alpha, A )

      call r8mat_print ( n, n, A, '  Gegenbauer to Monomial matrix G:' )  

      do i = 0, n - 1
        g(1:n) = 0.0D+00
        g(i+1) = 1.0D+00
        mono = matmul ( A, g )
        write ( *, '(a)' ) ''
        write ( *, '(a,i1)' ) 
     &    '  Monomial form of Gegenbauer polynomial ', i
        call r8poly_print ( i + 1, mono, '' )
      end do

      return
      end
      subroutine imtqlx_test ( )

c*********************************************************************72
c
cc imtqlx_test() tests imtqlx().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    10 June 2015
c
c  Author:
c
c    John Burkardt.
c
      implicit none

      integer, parameter :: n = 5

      double precision angle
      double precision d(n)
      double precision e(n)
      integer i
      double precision lam(n)
      double precision lam2(n)
      double precision qtz(n)
      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision z(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'imtqlx_test():'
      write ( *, '(a)' ) 
     &  '  imtqlx() takes a symmetric tridiagonal matrix A'
      write ( *, '(a)' ) '  and computes its eigenvalues LAM.'
      write ( *, '(a)' ) 
     &  '  It also accepts a vector Z and computes Q''*Z,'
      write ( *, '(a)' ) '  where Q is the matrix that diagonalizes A.'

      d(1:n) = 2.0D+00
      e(1:n-1) = -1.0D+00
      e(n) = 0.0D+00
      z(1:n) = 1.0D+00
c
c  On input, LAM is D, and QTZ is Z.
c
      lam(1:n) = d(1:n)
      qtz(1:n) = z(1:n)

      call imtqlx ( n, lam, e, qtz )

      call r8vec_print ( n, lam, '  Computed eigenvalues:' )

      do i = 1, n
        angle = dble ( i ) * r8_pi / dble ( 2 * ( n + 1 ) )
        lam2(i) = 4.0D+00 * ( sin ( angle ) ) ** 2
      end do

      call r8vec_print ( n, lam2, '  Exact eigenvalues:' )

      call r8vec_print ( n, z, '  Vector Z:' )
      call r8vec_print ( n, qtz, '  Vector Q''*Z:' )

      return
      end
      subroutine monomial_to_gegenbauer_matrix_test ( )

c*********************************************************************72
c
cc monomial_to_gegenbauer_matrix_test() tests monomial_to_gegenbauer_matrix().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    24 April 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 5

      double precision alpha
      double precision G(0:n-1,0:n-1)
      double precision Ginv(0:n-1,0:n-1)
      double precision I(0:n-1,0:n-1)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  monomial_to_gegenbauer_matrix_test():'
      write ( *, '(a)' ) 
     &  '  monomial_to_gegenbauer_matrix() evaluates the matrix'
      write ( *, '(a)' ) 
     &  '  which converts monomial polynomial coefficients'
      write ( *, '(a)' ) '  to Gegenbauer coefficients.'
      
      alpha = 0.5D+00

      call gegenbauer_to_monomial_matrix ( n, alpha, G )
      call r8mat_print ( n, n, G, '  Gegenbauer to Monomial matrix G:' )   

      call monomial_to_gegenbauer_matrix ( n, alpha, Ginv )
      call r8mat_print ( n, n, Ginv, 
     &  '  Monomial to Gegenbauer matrix Ginv:' )    

      I = matmul ( G, Ginv )
      call r8mat_print ( n, n, I, '  I = G * Ginv:' )    

      return
      end
      subroutine r8_hyper_2f1_test ( )

c*********************************************************************72
c
cc r8_hyper_2f1_test() tests r8_hyper_2f1().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    18 July 2015
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision a
      double precision b
      double precision c
      double precision fx
      double precision fx2
      integer n_data
      double precision x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'r8_hyper_2f1_test():'
      write ( *, '(a)' ) '  r8_hyper_2f1() evaluates the '
      write ( *, '(a)' ) '  hypergeometric 2F1 function.'
      write ( *, '(a)' ) ' '
      write ( *, '(a,a)' ) '      A       B       C       X      ', 
     &  ' 2F1                       2F1                     DIFF'
      write ( *, '(a,a)' ) '                                     ', 
     &  '(tabulated)               (computed)'
      write ( *, '(a)' ) ' '

      n_data = 0

      do

        call hyper_2f1_values ( n_data, a, b, c, x, fx )

        if ( n_data == 0 ) then
          exit
        end if

        call r8_hyper_2f1 ( a, b, c, x, fx2 )

        write ( *, 
     &    '(2x,f6.2,2x,f6.2,2x,f6.2,2x,f6.2,2x,' // 
     &    'g24.16,2x,g24.16,2x,g10.4)' ) 
     &    a, b, c, x, fx, fx2, abs ( fx - fx2 )

      end do

      return
      end
      subroutine r8_psi_test ( )

c*********************************************************************72
c
cc r8_psi_test() tests r8_psi().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    09 February 2008
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision fx
      double precision fx2
      integer n_data
      double precision r8_psi
      double precision x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'r8_psi_test():'
      write ( *, '(a)' ) '  r8_psi() evaluates the Psi function.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '      X         Psi(X)                     Psi(X)  ' //
     &  '               DIFF'
      write ( * , '(a)' ) 
     &  '               (Tabulated)                (R8_PSI)'
      write ( *, '(a)' ) ' '

      n_data = 0

      do

        call psi_values ( n_data, x, fx )

        if ( n_data == 0 ) then
          exit
        end if

        fx2 = r8_psi ( x )

        write ( *, '(2x,f8.4,2x,g24.16,2x,g24.16,2x,g10.4)' ) 
     &    x, fx, fx2, abs ( fx - fx2 )

      end do

      return
      end
      subroutine r8mat_print ( m, n, a, title )

c*********************************************************************72
c
cc r8mat_print() prints a real matrix.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 September 2004
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, the number of rows in A.
c
c    integer N, the number of columns in A.
c
c    double precision A(M,N), the matrix.
c
c    character ( len = * ) TITLE, a title.
c
      implicit none

      integer m
      integer n

      double precision a(m,n)
      character ( len = * ) title

      call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

      return
      end
      subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

c*********************************************************************72
c
cc r8mat_print_some prints some of a real matrix.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    10 September 2009
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, N, the number of rows and columns.
c
c    double precision A(M,N), an M by N matrix to be printed.
c
c    integer ILO, JLO, the first row and column to print.
c
c    integer IHI, JHI, the last row and column to print.
c
c    character ( len = * ) TITLE, a title.
c
      implicit none

      integer, parameter :: incx = 5
      integer m
      integer n

      double precision a(m,n)
      character ( len = 14 ) ctemp(incx)
      integer i
      integer i2hi
      integer i2lo
      integer ihi
      integer ilo
      integer inc
      integer j
      integer j2
      integer j2hi
      integer j2lo
      integer jhi
      integer jlo
      character ( len = * ) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )

      if ( m <= 0 .or. n <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  (None)'
        return
      end if

      do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

        j2hi = j2lo + incx - 1
        j2hi = min ( j2hi, n )
        j2hi = min ( j2hi, jhi )

        inc = j2hi + 1 - j2lo

        write ( *, '(a)' ) ' '

        do j = j2lo, j2hi
          j2 = j + 1 - j2lo
          write ( ctemp(j2), '(i8,6x)' ) j
        end do

        write ( *, '(''  Col   '',5a14)' ) ctemp(1:inc)
        write ( *, '(a)' ) '  Row'
        write ( *, '(a)' ) ' '

        i2lo = max ( ilo, 1 )
        i2hi = min ( ihi, m )

        do i = i2lo, i2hi

          do j2 = 1, inc

            j = j2lo - 1 + j2

            write ( ctemp(j2), '(g14.6)' ) a(i,j)

          end do

          write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc )

        end do

      end do

      return
      end
      subroutine r8poly_print ( n, a, title )

c*********************************************************************72
c
cc r8poly_print() prints a polynomial.
c
c  Discussion:
c
c    The power sum form is:
c
c      p(x) = a(0) + a(1) * x + ... + a(n-1) * x^(n-1) + a(n) * x^(n)
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    15 July 2015
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the dimension of A.
c
c    double precision A(0:N), the polynomial coefficients.
c    A(1) is the constant term and
c    A(N) is the coefficient of X^(N-1).
c
c    character ( len = * ) TITLE, a title.
c
      implicit none

      integer n

      double precision a(1:n)
      integer i
      double precision mag
      character plus_minus
      character ( len = * ) title

      if ( 0 < len_trim ( title ) ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) trim ( title )
      end if

      write ( *, '(a)' ) ' '

      if ( n < 0 ) then
        write ( *, '( ''  p(x) = 0'' )' )
        return
      end if

      if ( a(n) < 0.0D+00 ) then
        plus_minus = '-'
      else
        plus_minus = ' '
      end if

      mag = abs ( a(n) )

      if ( 3 <= n ) then
        write ( *, '( ''  p(x) = '', a1, g14.6, '' * x ^ '', i3 )' )
     &    plus_minus, mag, n - 1
      else if ( n == 2 ) then
        write ( *, '( ''  p(x) = '', a1, g14.6, '' * x'' )' ) 
     &    plus_minus, mag
      else if ( n == 1 ) then
        write ( *, '( ''  p(x) = '', a1, g14.6 )' ) plus_minus, mag
      end if

      do i = n - 1, 1, -1

        if ( a(i) < 0.0D+00 ) then
          plus_minus = '-'
        else
          plus_minus = '+'
        end if

        mag = abs ( a(i) )

        if ( mag /= 0.0D+00 ) then

          if ( 3 <= i ) then
            write ( *, 
     &        ' ( ''         '', a1, g14.6, '' * x ^ '', i3 )' ) 
     &        plus_minus, mag, i - 1
          else if ( i == 2 ) then
            write ( *, ' ( ''         '', a1, g14.6, '' * x'' )' ) 
     &        plus_minus, mag
          else if ( i == 1 ) then
            write ( *, ' ( ''         '', a1, g14.6 )' ) plus_minus, mag
          end if
        end if

      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the current YMDHMS date as a time stamp.
c
c  Example:
c
c    31 May 2001   9:45:54.872 AM
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    15 August 2021
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character ( len = 8 ) ampm
      integer d
      integer h
      integer m
      integer mm
      character ( len = 9 ), parameter, dimension(12) :: month = (/ 
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /)
      integer n
      integer s
      integer values(8)
      integer y

      call date_and_time ( values = values )

      y = values(1)
      m = values(2)
      d = values(3)
      h = values(5)
      n = values(6)
      s = values(7)
      mm = values(8)

      if ( h < 12 ) then
        ampm = 'AM'
      else if ( h == 12 ) then
        if ( n == 0 .and. s == 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h < 12 ) then
          ampm = 'PM'
        else if ( h == 12 ) then
          if ( n == 0 .and. s == 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end