!=======================================================================
! ch07_oscillatory_filon.f90
!
! Numerical integration of a rapidly oscillating function.
!
! The program evaluates
!
!          1
!          /
!         |  exp(-x) cos(omega*x) dx,       omega = 50,
!         /
!         0
!
! using two methods:
!
!   1. Composite Simpson rule
!   2. Composite quadratic Filon rule
!
! The Simpson rule resolves the complete oscillatory integrand
! numerically and therefore requires a relatively fine mesh.
!
! The Filon rule approximates only the slowly varying amplitude f(x)
! by a quadratic polynomial on pairs of subintervals.  The polynomial
! is then integrated analytically against cos(omega*x).  This allows a
! much coarser mesh when omega is large.
!
! This example follows Section 7.6, "Integration of rapidly oscillating
! functions," of Computational Physics: From Equations to Simulation.
!
! Alexander Godunov
! Companion website, 2026
!=======================================================================

program main
  implicit none

  double precision :: f
  double precision :: a, b, omega
  double precision :: integral, exact, error
  integer :: n, i
  integer, dimension(5) :: n_simpson
  integer, dimension(3) :: n_filon
  external :: f

  a     = 0.0d0
  b     = 1.0d0
  omega = 50.0d0

  ! Exact result for f(x) = exp(-x):
  !
  ! integral exp(-x) cos(omega*x) dx
  !   = exp(-x)[-cos(omega*x)+omega*sin(omega*x)]/(1+omega^2).
  !
  ! If f(x) is changed below, this reference value must also be
  ! replaced or omitted.
  exact = ( exp(-b)*(-cos(omega*b) + omega*sin(omega*b))        &
          - exp(-a)*(-cos(omega*a) + omega*sin(omega*a)) )      &
          /(1.0d0 + omega*omega)

  n_simpson = (/ 100, 200, 400, 800, 1600 /)
  n_filon   = (/  20,  40,  80 /)

  write(*,'(a)')        'Rapidly oscillating integral'
  write(*,'(a,f6.1)')   'Integral of exp(-x) cos(omega*x), omega = ', omega
  write(*,'(a,es15.7)') 'Exact value          = ', exact

  write(*,*)
  write(*,'(a)') 'Composite Simpson rule'
  write(*,'(a)') '        n            Integral       Absolute error'
  write(*,'(a)') '----------------------------------------------------'

  do i = 1, size(n_simpson)

    n = n_simpson(i)
    call composite_simpson_oscillatory(f, a, b, omega, n, integral)

    error = abs(integral-exact)
    write(*,'(i9,2es20.7)') n, integral, error

  end do

  write(*,*)
  write(*,'(a)') 'Composite quadratic Filon rule'
  write(*,'(a)') '        n            Integral       Absolute error'
  write(*,'(a)') '----------------------------------------------------'

  do i = 1, size(n_filon)

    n = n_filon(i)
    call filon_quadratic(f, a, b, omega, n, integral)

    error = abs(integral-exact)
    write(*,'(i9,2es20.7)') n, integral, error

  end do

end program main


double precision function f(x)
  !---------------------------------------------------------------------
  ! Slowly varying amplitude in the oscillatory integral
  !
  !          b
  !          /
  !         |  f(x) cos(omega*x) dx .
  !         /
  !          a
  !
  ! To use this program for another problem with the same oscillatory
  ! factor, replace only the line defining f(x) below.  The Simpson and
  ! Filon routines do not need to be changed.
  !
  ! If an analytic reference value is used in the main program, it must
  ! also be changed when f(x) is changed.
  !
  ! Input:
  !   x - integration variable
  !
  ! Output:
  !   f - value of the slowly varying amplitude
  !---------------------------------------------------------------------
  implicit none

  double precision, intent(in) :: x

  f = exp(-x)

end function f


subroutine composite_simpson_oscillatory(f, a, b, omega, n, integral)
  !---------------------------------------------------------------------
  ! Composite Simpson rule for
  !
  !          b
  !          /
  !         |  f(x) cos(omega*x) dx .
  !         /
  !          a
  !
  ! Input:
  !   f     - slowly varying amplitude
  !   a     - lower integration limit
  !   b     - upper integration limit
  !   omega - angular frequency of the oscillatory factor
  !   n     - number of equal subintervals; n must be even
  !
  ! Output:
  !   integral - numerical approximation to the integral
  !
  ! The complete oscillatory integrand is sampled numerically.  For
  ! large omega, n must therefore be large enough to resolve the rapid
  ! oscillations.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b, omega
  integer, intent(in)           :: n
  double precision, intent(out) :: integral
  external :: f

  double precision :: h, x, sum_odd, sum_even
  integer :: i

  if (n < 2 .or. mod(n,2) /= 0) then
    write(*,'(a)') 'composite_simpson_oscillatory: n must be even.'
    stop
  end if

  h = (b-a)/dble(n)

  sum_odd = 0.0d0
  do i = 1, n-1, 2
    x = a + dble(i)*h
    sum_odd = sum_odd + f(x)*cos(omega*x)
  end do

  sum_even = 0.0d0
  do i = 2, n-2, 2
    x = a + dble(i)*h
    sum_even = sum_even + f(x)*cos(omega*x)
  end do

  integral = h*( f(a)*cos(omega*a) + f(b)*cos(omega*b)          &
             + 4.0d0*sum_odd + 2.0d0*sum_even )/3.0d0

end subroutine composite_simpson_oscillatory


subroutine filon_quadratic(f, a, b, omega, n, integral)
  !---------------------------------------------------------------------
  ! Composite quadratic Filon rule for
  !
  !          b
  !          /
  !         |  f(x) cos(omega*x) dx .
  !         /
  !          a
  !
  ! Input:
  !   f     - slowly varying amplitude
  !   a     - lower integration limit
  !   b     - upper integration limit
  !   omega - angular frequency of the oscillatory factor
  !   n     - number of equal subintervals; n must be even
  !
  ! Output:
  !   integral - numerical approximation to the integral
  !
  ! Method:
  !   The subintervals are taken in pairs.  On each pair
  !
  !       [x0,x2],       x1 = (x0+x2)/2,
  !
  !   introduce the local coordinate
  !
  !       x = x1 + t,       -h <= t <= h.
  !
  !   The slowly varying amplitude is approximated by
  !
  !       f(x1+t) = a0 + a1*t + a2*t^2,
  !
  !   where the coefficients are obtained from f(x0), f(x1), f(x2).
  !
  !   The oscillatory factor is NOT approximated.  The polynomial is
  !   integrated analytically against cos[omega(x1+t)].
  !
  !   By symmetry over [-h,h],
  !
  !       integral cos(omega*t) dt        = C0,
  !       integral t^2 cos(omega*t) dt    = C2,
  !       integral t sin(omega*t) dt      = S1,
  !
  !   while
  !
  !       integral t cos(omega*t) dt      = 0,
  !       integral sin(omega*t) dt        = 0,
  !       integral t^2 sin(omega*t) dt    = 0.
  !
  !   Thus the contribution from one pair is
  !
  !     cos(omega*x1)*(a0*C0 + a2*C2)
  !       - sin(omega*x1)*a1*S1 .
  !
  !   This is the essential Filon idea: approximate the smooth
  !   amplitude, but integrate the rapid oscillation analytically.
  !---------------------------------------------------------------------
  implicit none

  double precision :: f
  double precision, intent(in)  :: a, b, omega
  integer, intent(in)           :: n
  double precision, intent(out) :: integral
  external :: f

  double precision :: h
  double precision :: x0, x1, x2
  double precision :: f0, f1, f2
  double precision :: a0, a1, a2
  double precision :: theta
  double precision :: c0, c2, s1
  double precision :: pair_integral
  integer :: k

  if (n < 2 .or. mod(n,2) /= 0) then
    write(*,'(a)') 'filon_quadratic: n must be a positive even integer.'
    stop
  end if

  if (abs(omega) < 1.0d-12) then
    write(*,'(a)') 'filon_quadratic: omega is too close to zero.'
    stop
  end if

  h = (b-a)/dble(n)
  theta = omega*h

  ! Analytic moments over -h <= t <= h.
  c0 = 2.0d0*sin(theta)/omega

  c2 = 2.0d0*h*h*sin(theta)/omega                         &
     + 4.0d0*h*cos(theta)/(omega*omega)                   &
     - 4.0d0*sin(theta)/(omega*omega*omega)

  s1 = -2.0d0*h*cos(theta)/omega                          &
     +  2.0d0*sin(theta)/(omega*omega)

  integral = 0.0d0

  do k = 0, n-2, 2

    x0 = a + dble(k)*h
    x1 = x0 + h
    x2 = x0 + 2.0d0*h

    f0 = f(x0)
    f1 = f(x1)
    f2 = f(x2)

    ! Quadratic interpolation of the amplitude in local coordinate t.
    a0 = f1
    a1 = (f2-f0)/(2.0d0*h)
    a2 = (f0-2.0d0*f1+f2)/(2.0d0*h*h)

    pair_integral = cos(omega*x1)*(a0*c0 + a2*c2)         &
                  - sin(omega*x1)*a1*s1

    integral = integral + pair_integral

  end do

end subroutine filon_quadratic