[nichollsh]

Description

SIMPSONS_RULE assumed an even number of intervals. With an odd number, the last interval counted at about two thirds of its value. CALC_PLANCK_TBL passes an odd number for about half of all bands and temperatures. As a result, the Planck table that prep_spec writes in block 6 was too low. For a file with 4096 bands of equal width, the error was 0.5 % at 2604 K and 0.9 % at 4494 K.

With an odd number of intervals, the routine now uses Simpson's 3/8 rule over the last three intervals. A single interval uses the trapezoidal rule. Results for an even number of intervals do not change. One side effect: N = 0 now returns 0 instead of INPUT_ARRAY(0)*DELTA/3.

The filter function and polynomial fit options of block 6 call the same routine, so they also get the fix. Ccorr_k does not use this routine. The version is now 2607.3.

Closes #35

Validation of changes

Test configuration: Fedora 44, x86-64, GNU Fortran 16.2.1. Built with ./configure && ./build_code (-Ofast -march=native). The build has no warnings for simpsons_rule.f.

Integrator alone: I compiled a small driver against the old and new simpsons_rule.f and integrated over [0, 1] with N = 1 to 12 intervals.

• Constant function: for odd N the old result is low by 1/(3N), for example -0.111 at N = 3 and -0.030 at N = 11. The new result is exact for every N.
• x^3: the new result is exact for every N >= 2.
• exp(x): the new error at odd N is in line with the neighbouring even N, for example 1.5e-4 at N = 3 and 5.6e-7 at N = 11.
• Every even N gives the same result as before.

Planck table: I ran the old and new prep_spec on copies of three spectral files, with the same inputs that AGNI uses:

printf '%s\n' file.sf a 6 n T "25 6000" 1196 -1 | $RAD_DIR/bin/prep_spec

I then compared each band of the table with 64-point Gauss-Legendre quadrature of the Planck function. 128-point quadrature agrees to 1e-14.

Sum of the table over all bands, divided by sigma T^4, for a file with 4096 bands from 10 to 50000 cm-1:

T [K]	Old	New	Emission inside the file's range
998.5	0.99809	0.99999	1.00000
1796	0.99655	0.99999	1.00000
2604	0.99502	0.99999	1.00000
3002	0.99999	0.99999	1.00000
4494	0.99136	0.99990	0.99991

Largest relative difference from quadrature, over bands that carry more than 1e-6 of sigma T^4, at the same five temperatures:

File	Old	New
4096 bands, 10-50000 cm-1	up to 8.9e-3	up to 9.0e-5
64 bands, 10-50000 cm-1	up to 1.9e-4	up to 1.2e-4
Honeyside, 48 bands	up to 3.4e-3	up to 1.2e-4

The remaining difference of about 1e-4 also appears in the old build where the interval count was already even, so the integration does not cause it. I did not trace it.

Regression tests: I started examples/quick_tests with this build but it did not finish. The first test (runes) reported Files output.txt and gfortran_8_1_0.txt differ, and the same for the ifort reference. I have not checked whether the build without this change gives the same differences. A later test then looped on +++ Unexpected response:, and I stopped the run. I did not run examples/slow_tests, the Fortitude lint, or a full AGNI run with this build.

Checklist

• I have followed the contributing guidelines
• My code follows the style guidelines of this project
• I have performed a self-review of my code
• My changes generate no new warnings or errors
• I have checked that the tests still pass on my computer
• I have updated the docs, as appropriate
• I have added tests for these changes, as appropriate
• I have checked that all dependencies have been updated, as required