//=======================================================================
// ch07_oscillatory.cpp
//
// 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
//
// Simpson resolves the complete oscillatory integrand numerically.
// Filon approximates only the slowly varying amplitude by a quadratic
// polynomial and integrates the oscillatory factor analytically.
//
// Based on the Fortran implementation by Alexander Godunov.
// C++ version prepared for the companion website, 2026.
//=======================================================================

#include <cmath>
#include <iomanip>
#include <iostream>

using namespace std;

double f(double x);

double composite_simpson_oscillatory(double (*fun)(double),
                                     double a, double b,
                                     double omega, int n);

double filon_quadratic(double (*fun)(double),
                       double a, double b,
                       double omega, int n);

int main()
{
    const double a = 0.0;
    const double b = 1.0;
    const double omega = 50.0;

    // Exact result for f(x) = exp(-x).
    const double exact =
        (exp(-b)*(-cos(omega*b) + omega*sin(omega*b))
       - exp(-a)*(-cos(omega*a) + omega*sin(omega*a)))
       /(1.0 + omega*omega);

    const int n_simpson[] = {100, 200, 400, 800, 1600};
    const int n_filon[]   = {20, 40, 80};

    cout << "Rapidly oscillating integral\n";
    cout << "Integral of exp(-x) cos(omega*x), omega = "
         << fixed << setprecision(1) << omega << '\n';
    cout << scientific << setprecision(7);
    cout << "Exact value          = " << setw(15) << exact << "\n\n";

    cout << "Composite Simpson rule\n";
    cout << "        n            Integral       Absolute error\n";
    cout << "----------------------------------------------------\n";

    for (int n : n_simpson) {
        const double integral =
            composite_simpson_oscillatory(f, a, b, omega, n);
        const double error = fabs(integral-exact);

        cout << setw(9) << defaultfloat << n
             << scientific << setprecision(7)
             << setw(20) << integral
             << setw(20) << error << '\n';
    }

    cout << "\nComposite quadratic Filon rule\n";
    cout << "        n            Integral       Absolute error\n";
    cout << "----------------------------------------------------\n";

    for (int n : n_filon) {
        const double integral =
            filon_quadratic(f, a, b, omega, n);
        const double error = fabs(integral-exact);

        cout << setw(9) << defaultfloat << n
             << scientific << setprecision(7)
             << setw(20) << integral
             << setw(20) << error << '\n';
    }

    return 0;
}


double f(double x)
{
    //-------------------------------------------------------------------
    // Slowly varying amplitude in
    //
    //       integral f(x) cos(omega*x) dx.
    //
    // To use another amplitude, replace only this line.
    //-------------------------------------------------------------------
    return exp(-x);
}


double composite_simpson_oscillatory(double (*fun)(double),
                                     double a, double b,
                                     double omega, int n)
{
    //-------------------------------------------------------------------
    // Composite Simpson rule for f(x) cos(omega*x).
    //-------------------------------------------------------------------

    if (n < 2 || n % 2 != 0) {
        cerr << "composite_simpson_oscillatory: n must be even.\n";
        return 0.0;
    }

    const double h = (b-a)/static_cast<double>(n);

    double sum_odd = 0.0;
    for (int i = 1; i < n; i += 2) {
        const double x = a + static_cast<double>(i)*h;
        sum_odd += fun(x)*cos(omega*x);
    }

    double sum_even = 0.0;
    for (int i = 2; i <= n-2; i += 2) {
        const double x = a + static_cast<double>(i)*h;
        sum_even += fun(x)*cos(omega*x);
    }

    return h*(fun(a)*cos(omega*a)
            + fun(b)*cos(omega*b)
            + 4.0*sum_odd + 2.0*sum_even)/3.0;
}


double filon_quadratic(double (*fun)(double),
                       double a, double b,
                       double omega, int n)
{
    //-------------------------------------------------------------------
    // Composite quadratic Filon rule for
    //
    //       integral f(x) cos(omega*x) dx.
    //
    // On each pair of subintervals, f(x) is approximated by a
    // quadratic polynomial in a local coordinate. The polynomial is
    // then integrated analytically against the oscillatory factor.
    //-------------------------------------------------------------------

    if (n < 2 || n % 2 != 0) {
        cerr << "filon_quadratic: n must be a positive even integer.\n";
        return 0.0;
    }

    if (fabs(omega) < 1.0e-12) {
        cerr << "filon_quadratic: omega is too close to zero.\n";
        return 0.0;
    }

    const double h = (b-a)/static_cast<double>(n);
    const double theta = omega*h;

    // Analytic moments over -h <= t <= h.
    const double c0 = 2.0*sin(theta)/omega;

    const double c2 =
          2.0*h*h*sin(theta)/omega
        + 4.0*h*cos(theta)/(omega*omega)
        - 4.0*sin(theta)/(omega*omega*omega);

    const double s1 =
         -2.0*h*cos(theta)/omega
        + 2.0*sin(theta)/(omega*omega);

    double integral = 0.0;

    for (int k = 0; k <= n-2; k += 2) {

        const double x0 = a + static_cast<double>(k)*h;
        const double x1 = x0 + h;
        const double x2 = x0 + 2.0*h;

        const double f0 = fun(x0);
        const double f1 = fun(x1);
        const double f2 = fun(x2);

        const double a0 = f1;
        const double a1 = (f2-f0)/(2.0*h);
        const double a2 = (f0-2.0*f1+f2)/(2.0*h*h);

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

        integral += pair_integral;
    }

    return integral;
}