◆ deoint_r()

Semi-infinite interval automatic quadrature for Fourier integrals (double exponential (DE) formula) (reverse communication version)

Purpose: This routine computes the following Fourier integrals by the automatic integration applying the double exponential (DE) formula.

Ic = ∫ f(x) cosωx dx [a, +∞] or Is = ∫ f(x) sinωx dx [a, +∞]

where f(x) is computed and provided by the user in accordance with irev.

Ic or Is is transformed to the intergral on the infinite interval [-∞, +∞] by using the following transformation function, and the integral value is computed by the trapezoidal rule.
Ic: x = Mφ(t - π/2M)/ω

Is: x = Mφ(t)/ω

The following φ(t) is used.
φ1(t) = t/(1 - exp(-2t -α(1 - exp(-t)) - β(exp(t) - 1))),

β = 1/4, α = β/sqrt(1 + Mlog(1 + M)/4π)

or

φ2(t) = t/(1 - exp(-ksinh(t))), k = 2π

The step size h and the constant M are selected to satisfy Mh = π.
Please refer to the reference below for more details.

Parameters

[in]	a	Lower bound of integration interval.
[in]	omega	ω value in cos or sin.
[in]	iflag	Choose the type of integral. = 0: ∫ f(x) cosωx dx [0, +∞] = 1: ∫ f(x) sinωx dx [0, +∞] Further, the transformation function can be selected as follows. iflag += 0: φ1(t) (recommended) iflag += 2: φ2(t)
[in]	eps	Relative accuracy requested. (0 < eps <= 1)
[out]	result	Obtained approximation to the integral.
[out]	neval	Number of evaluations of f(x).
[out]	err	Estimated relative error of obtained v. To be overestimated for many cases. This routine may return as successful exit even if err > eps.
[out]	info	= 0: Successful exit = -3: The argument iflag had an illegal value (iflag != 0 to 3) = -4: The argument eps had an illegal value (eps <= 0 or eps > 1) = 1: Not converged
[out]	xx	irev = 1, 2: xx contains the abscissa where the function value should be evaluated and given in the next call.
[in]	yy	irev = 1, 2: The function value f(xx) should be given in yy in the next call.
[in,out]	irev	Control variable for reverse communication. [in] Before first call, irev should be initialized to zero. On succeeding calls, irev should not be altered. [out] If irev is not zero, complete the following tasks and call this routine again without changing irev. = 0: Computation finished. = 1, 2: User should set the function value at xx in yy. Do not alter any variables other than yy.

Reference

Ooura and Mori, "The double exponential formula for oscillatory functions over the half infinite interval", J. Comput. Appl. Math., 38 (1991), 353-360.
Ooura and Mori, "A robust double exponential formula for Fourier-type integrals", J. Comput. Appl. Math., 112 (1999), 229-241.