================================================================================
REPRODUCIBILITY PACKAGE — Journal of Engineering Mathematics submission
"Perturbative modelling of Arrhenius reactive transport in non-isothermal
 Gordon-Schowalter Couette-Poiseuille flow"
Soria and Harrup — v2-clean, June 2026
================================================================================

REQUIREMENTS
  Python >= 3.10 with numpy, scipy and matplotlib.

SCRIPTS
  Run each script from this directory. No command-line arguments are required.

1) validate.py
   Full numerical validation of the manuscript:
     S0  Sign adjudication of the viscoelastic thermal correction using the
         exact rational Johnson-Segalman/Gordon-Schowalter shear law.
     S1  Sturm-Liouville shooting solver cross-check against an independent
         finite-difference/generalised-eigenvalue implementation.
     S2  Closed-form temperature fields theta_B and theta_BW vs direct
         double-quadrature integration.
     S3  Cross-term validation: exact perturbed eigenvalue vs first-order
         formula (eqs. 14-15) over the (a, Ws) grid of Table 2; route
         additivity; UCM/LCM invariance; crossover number K.
     S3b Crossover criterion sweep corresponding to Table 3.
   Outputs: report.txt, results.json.

2) fig_gen.py
   Regenerates Figures 1-3 and the data behind Table 2; includes an independent
   second-order finite-difference cross-check of the shooting solver
   (generalised eigenproblem, 4001 nodes).
   Outputs: fig1_schematic.png, fig2_thermal.png, fig3_collapse.png,
   figdata.json.

3) fig4_table3.py
   Computes the K design map on the thermally corrected Newtonian base state,
   draws the Frank-Kamenetskii validity boundary gamma_A*Theta_max = 1, and
   reproduces/extends Table 3.
   Outputs: fig4_Kmap.png, table3.json.

TOLERANCES
  ODE integration: adaptive RK45, rtol 1e-11, atol 1e-13.
  Eigenvalue refinement: Brent, xtol 1e-12.
  FD cross-check grid: N = 4001 nodes.

KEY NUMBERS (operating point Pe_M=5, Da_II=1, P*=3, Br=1, gamma_A=5)
  lambda_1 = 1.28025  (lambda_1^2 = 1.63905)
  Newtonian thermal shift  dl_B = +3.391e-2
  Route coefficients per unit eps: thermal-reactive +3.124e-2,
                                   advective        -5.003e-1
  K = 0.0625; worst first-order error about 2.2%; route additivity about 0.13%;
  UCM/LCM shift = 0 to machine precision.

AUTHORSHIP NOTE
  These files correspond to the two-author version of the manuscript:
  Leonardo D. Soria R. and Anthony A. Harrup G.
================================================================================