#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
fig4_table3.py
===================
Crossover-map and Table 3 script for the manuscript:

  "Perturbative modelling of Arrhenius reactive transport in non-isothermal
   Gordon-Schowalter Couette-Poiseuille flow"

Generates:
  - Figure 4: design map of the crossover number K.
  - Frank-Kamenetskii validity boundary gamma_A*Theta_max = 1.
  - table3.json: crossover-verification data.

Authors: Leonardo D. Soria R. and Anthony A. Harrup G.
Target journal: Journal of Engineering Mathematics
Version: v2-clean, June 2026
"""
import numpy as np, json
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
plt.rcParams.update({"font.family":"serif","font.size":10,"figure.dpi":200})
RNG=dict(method="RK45",rtol=1e-11,atol=1e-13)
P,Pe=3.0,5.0
phi0=lambda e:1.0+(P/2)*e*(1-e)
phi1=lambda e:(P**3/4)*(1/16-(e-0.5)**4)
thB =lambda e:(P**2/12)*(1/16-(e-0.5)**4)
thBW=lambda e:(P**4/30)*(1/64-(e-0.5)**6)
def res(lam,Da_f,phi_f):
    r=lambda t,y:[y[1],(Da_f(t)-lam**2*Pe*phi_f(t))*y[0]]
    return solve_ivp(r,[0.5,1.0],[1.0,0.0],**RNG).y[0,-1]
def lam1(Da_f,phi_f,lmax=10,dl=0.05):
    g=np.arange(dl,lmax,dl); rp=res(g[0],Da_f,phi_f)
    for l in g[1:]:
        r=res(l,Da_f,phi_f)
        if rp*r<0: return brentq(res,l-dl,l,args=(Da_f,phi_f),xtol=1e-12,rtol=1e-13)
        rp=r
    raise RuntimeError
def mode(lam,Da_f,phi_f,N=2001):
    eta=np.linspace(0,1,N); half=eta[eta>=0.5]
    r=lambda t,y:[y[1],(Da_f(t)-lam**2*Pe*phi_f(t))*y[0]]
    s=solve_ivp(r,[0.5,1.0],[1.0,0.0],t_eval=half,**RNG)
    H=np.concatenate([s.y[0][:0:-1],s.y[0]]); return eta,H/np.max(np.abs(H))
def Kcorr(Da0,BG):
    Da_b=lambda e: Da0*(1.0+BG*thB(e))
    lb=lam1(Da_b,phi0); eta,Hb=mode(lb,Da_b,phi0)
    I=lambda f: np.trapezoid(f,eta)
    return BG*Da0*I(thBW(eta)*Hb**2)/(lb**2*Pe*I(phi1(eta)*Hb**2)), lb
def net_shift(Da0,BG,eps=0.0625):
    Da_b=lambda e: Da0*(1.0+BG*thB(e))
    Da_p=lambda e: Da0*(1.0+BG*(thB(e)+eps*thBW(e)))
    phi_p=lambda e: phi0(e)+eps*phi1(e)
    return lam1(Da_p,phi_p)**2-lam1(Da_b,phi0)**2

# --- Table 3 (reproduce v1 + new in-validity rows) ---
print("Tabla 3 (a=0, Ws=0.25, eps=0.0625, Pe=5, P*=3):")
rows=[]
for Da0,BG in [(1,2),(1,60),(5,30),(5,60),(10,12),(10,20)]:
    K,_=Kcorr(Da0,BG); dl=net_shift(Da0,BG)
    rows.append(dict(Da=Da0,BrgA=BG,K=K,dl=dl,gATmax=BG*P**2/192))
    print(f"  Da={Da0:>2} BrgA={BG:>3}  K={K:.3f}  dl_net={dl:+.3e}  sign={'+' if dl>0 else '-'}  gA*Tmax={BG*P**2/192:.2f}")

# --- K map with corrected base + validity line ---
Das=np.logspace(np.log10(0.3),np.log10(12),13)
BGs=np.logspace(np.log10(0.5),np.log10(100),11)
KK=np.zeros((len(BGs),len(Das)))
for i,bg in enumerate(BGs):
    for j,da in enumerate(Das):
        KK[i,j]=Kcorr(da,bg)[0]
DD,GG=np.meshgrid(Das,BGs)
fig,ax=plt.subplots(figsize=(4.8,3.5))
cs=ax.contourf(DD,GG,np.log10(KK),levels=np.linspace(-2.5,0.8,14),cmap="RdBu_r")
c1=ax.contour(DD,GG,KK,levels=[1.0],colors="k",linewidths=2.0)
ax.clabel(c1,fmt={1.0:"K = 1"},fontsize=9)
c2=ax.contour(DD,GG,KK,levels=[0.1],colors="k",linewidths=0.9,linestyles="--")
ax.clabel(c2,fmt={0.1:"K = 0.1"},fontsize=8)
ax.axhline(192/P**2,color="0.25",lw=1.2,ls=":")
ax.text(0.33,192/P**2*1.12,r"$\gamma_A\Theta_{\max}=1$ (FK linearisation boundary)",fontsize=7.5)
for r in rows:
    ax.plot(r["Da"],r["BrgA"],"ks" if r["dl"]>0 else "ko",ms=7,
            mfc=("k" if r["dl"]>0 else "none"))
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel(r"$Da_{II}$"); ax.set_ylabel(r"$Br\,\gamma_A$")
ax.set_title(r"Design map of $K$ ($Pe_M=5$, $P^*=3$)",fontsize=10)
cb=plt.colorbar(cs,ax=ax); cb.set_label(r"$\log_{10}K$",fontsize=9)
ax.text(0.35,1.1,"advection-dominated\n($K<1$): longer reactor",fontsize=7.5)
ax.text(4.0,55,"kinetics-\ndominated\n($K>1$)",fontsize=7.5,ha="center")
plt.tight_layout(); plt.savefig("fig4_Kmap.png",bbox_inches="tight"); plt.close()
json.dump(dict(table3=rows),open("table3.json","w"),indent=1)
print("fig4 regenerada (base corregida) + table3.json")