subroutine sine_transform_data ( n, d, s )

!*****************************************************************************80
!
!! sine_transform_data() does a sine transform on a vector of data.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 February 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data points.
!
!    Input, real ( kind = rk ) D(N), the vector of data.
!
!    Output, real ( kind = rk ) S(N), the sine transform coefficients.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) angle
  real ( kind = rk ) d(n)
  integer i
  integer j
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) s(n)

  do i = 1, n
    s(i) = 0.0D+00
    do j = 1, n
      angle = pi * real ( i * j, kind = rk ) / real ( n + 1, kind = rk )
      s(i) = s(i) + sin ( angle ) * d(j)
    end do
    s(i) = s(i) * sqrt ( 2.0D+00 / real ( n + 1, kind = rk ) )
  end do

  return
end
subroutine sine_transform_function ( n, a, b, f, s )

!*****************************************************************************80
!
!! sine_transform_function() does a sine transform on functional data.
!
!  Discussion:
!
!    The interval [A,B] is divided into N+1 intervals using N+2 points,
!    which are indexed by 0 through N+1.
!
!    The original function F(X) is regarded as the sum of a linear function 
!    F1 that passes through (A,F(A)) and (B,F(B)), and a function F2
!    which is 0 at A and B.
!
!    The sine transform coefficients for F2 are then computed.
!
!    To recover the interpolant of F(X), it is necessary to combine the
!    linear part F1 with the sine transform interpolant:
!
!      Interp(F)(X) = F1(X) + F2(X)
!
!    This can be done by calling SINE_TRANSFORM_INTERPOLANT().
!    
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 February 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data points.
!
!    Input, real ( kind = rk ) A, B, the interval endpoints.
!
!    Input, external, real ( kind = rk ) F, a pointer to the function.
!
!    Output, real ( kind = rk ) S(N), the sine transform coefficients.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a
  real ( kind = rk ) angle
  real ( kind = rk ) b
  real ( kind = rk ), external :: f
  real ( kind = rk ) f2(n)
  real ( kind = rk ) fa
  real ( kind = rk ) fb
  integer i
  integer j
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) s(n)
  real ( kind = rk ) x(n)
!
!  Evenly spaced points between A and B, but omitting
!  A and B themselves.
!
  do i = 1, n
    x(i) = ( real ( n - i + 1, kind = rk ) * a   &
           + real (     i,     kind = rk ) * b ) &
           / real ( n     + 1, kind = rk )
  end do
!
!  Subtract F1(X) from F(X) to get F2(X).
!
  fa = f ( a )
  fb = f ( b )

  do i = 1, n
    f2(i) = f ( x(i) )                &
          - ( ( b - x(i)     ) * fa   &
          +   (     x(i) - a ) * fb ) &
          /   ( b        - a )
  end do
!
!  Compute the sine transform of F2(X).
!
  do i = 1, n
    s(i) = 0.0D+00
    do j = 1, n
      angle = pi * real ( i * j, kind = rk ) / real ( n + 1, kind = rk )
      s(i) = s(i) + sin ( angle ) * f2(j)
    end do
    s(i) = s(i) * sqrt ( 2.0D+00 / real ( n + 1, kind = rk ) )
  end do

  return
end
subroutine sine_transform_interpolant ( n, a, b, fa, fb, s, nx, x, value )

!*****************************************************************************80
!
!! sine_transform_interpolant() evaluates the sine transform interpolant.
!
!  Discussion:
!
!    The interval [A,B] is divided into N+1 intervals using N+2 points,
!    which are indexed by 0 through N+1.
!
!    The original function F(X) is regarded as the sum of a linear function 
!    F1 that passes through (A,F(A)) and (B,F(B)), and a function F2
!    which is 0 at A and B.
!
!    The function F2 has been approximated using the sine transform,
!    and the interpolant is then evaluated as:
!
!      Interp(F)(X) = F1(X) + F2(X)
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 February 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of terms in the approximation.
!
!    Input, real ( kind = rk ) A, B, the interval over which the approximant 
!    was defined.
!
!    Input, real ( kind = rk ) FA, FB, the function values at A and B.
!
!    Input, real ( kind = rk ) S(N), the approximant coefficients.
!
!    Input, integer NX, the number of evaluation points.
!
!    Input, real ( kind = rk ) X(NX), the evaluation points.
!
!    Output, real ( kind = rk ) VALUE(NX), the value of the interpolant.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  integer nx

  real ( kind = rk ) a
  real ( kind = rk ) angle
  real ( kind = rk ) b
  real ( kind = rk ) f1(nx)
  real ( kind = rk ) f2(nx)
  real ( kind = rk ) fa
  real ( kind = rk ) fb
  integer i
  integer j
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) s(n)
  real ( kind = rk ) value(nx)
  real ( kind = rk ) x(nx)
!
!  Compute linear function F1(X).
!
  f1(1:nx) = ( ( b - x(1:nx)     ) * fa   &
             + (     x(1:nx) - a ) * fb ) &
             / ( b           - a )
!
!  Compute sine interpolant F2(X).
!
  f2(1:nx) = 0.0D+00

  do i = 1, nx
    do j = 1, n
      angle = real ( j, kind = rk ) * ( x(i) - a ) * pi / ( b - a )
      f2(i) = f2(i) + s(j) * sin ( angle )
    end do
  end do

  f2(1:nx) = f2(1:nx) * sqrt ( 2.0D+00 / real ( n + 1, kind = rk ) )
!
!  Interpolant = F1 + F2.
! 
  value(1:nx) = f1(1:nx) + f2(1:nx)

  return
end