Caution When Using Adaptive Integrators
			M.Lampton UCB SSL 2009

In optics, a common task is to evaluate a system's modulation transfer
function (MTF) from component factors, and then Bessel transform that MTF
to evaluate the point-spread function (PSF) and encircled energy (EE). 

In this work I usually use an adaptive integrator: adaptive here in the
sense that rather than simply pre-specifying (say) 1000 equally spaced
samples on the (a,b) integration axis, and a fixed trapezoidal rule, 
I instead use a routine like Press et al ("Numerical Recipes") QSIMP
that starts out with a coarse sampling and iterates with successively
finer samples until the results stabilize to a pre-specified tolerance.

I have found that sometimes the MTF-to-PSF integrator fails to "catch" 
--- it exits early, returning zero rather than its correct PSF value.

This is particularly apparent with short wavelengths and large damping
attenuation of the MTF for example with large added Gaussian or exponential
blur dialed in, which are low Strehl circumstances.

This premature exit has two causes:

(1) The integrand is periodic and rapidly decaying.  The rapid decay
comes from the low Strehl, so the integrand is concentrated near v=0.
The periodic part is of course the J0 Bessel function, similar to a cosine.

(2) The succession of binary subdivisions of the (a,b) axis is permitted
to exit with sparse sampling of the (a,b) interval.  So if only 128
samples are chosen, and if the first few of these are all near the zeroes
of the Bessel function, then the routine will exit early with a near-zero 
incorrect result. If the integrand lies concentrated near the lower limit,
and if it has a few zeros there, there will be a band of false zero exits. 

This situation has several possible remedies:

(1) Simply jack up the MINITER threshold so that sufficient sampling will always
take place.  For example with MINITER=8, doRomberg() takes 257 samples including
both end points.  Or with MINITER=12, takes 4097 samples minimum.   If your
integrand has only 100 zeros in (a,b) then MINITER=8 should be plenty.

(2) Or tighten up the exit tolerance TOL, making it (say) 1E-12 or 1E-14, 
so that even a hint of an imperfectly zero sample will force the iterations
to continue.  The downside is that excessive accuracy demands more machine
time & iterations, even when no sampling resonances are present. 

(3) Or more radically, adopt a nonperiodic sample method.  If you know
in advance that most of the integrand lies near "a" and fades rapidly towards 
"b" then equal density samping of the (a,b) interval is computationally wasteful.  
I offer an example of a nonuniform integrator "Nonuniform Integration" on my website.
There, I map a uniform sampling interval 0<u<1 onto the integration variable
0<x<1 with a parabolic relation, x=u^2.  This concentrates the bottom 10% of the
samples onto the bottom 1% of the integration range, which is a big help
in establishing accurate convergence when the integrand is peaked at x=a.
Moreover, the periodic samples in u are nonperiodic in x, so that an integrand
that is oscillatory in x cannot resonate with its samples.