"""Engineering judgment cases, selected without fitting historical residues."""
import math,json,hashlib
from pathlib import Path
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
q=math.exp(8205.7/284.15-25.323)*.02*.03*26.018*1000
ref=101.5/8800*10000
tau=np.array([7.2,96.]);weights=100/(-np.expm1(-24/tau));weights/=sum(weights)
g=10/(.02*101325*27.026/(8.31446261815324*284.15))
def integrate(k,n,duration,cycles,level):
 y=np.zeros(4)
 def step(y,h,c):
  def f(t,y):
   m=np.maximum(y[:2],0);total=sum(m);r=k*total**n
   capture=r*m/total if total>0 else np.zeros(2)
   return np.r_[q*weights*c/tau-m/tau-capture,r,sum(m/tau)]
  return solve_ivp(f,[0,h],y,method='DOP853',rtol=1e-9,atol=1e-10).y[:,-1]
 for i in range(cycles):
  y=step(y,duration,level)
  if i+1<cycles:y=step(y,24-duration,0)
 y=step(y,1704,0)
 incoming=q*level*duration*cycles*sum(weights/tau)
 error=abs(sum(y)-incoming)
 assert error<1e-5 and min(y)>-1e-7
 assert y[2]<10000*18*26.018/(7*55.845)
 return dict(retained_CN_mg_kg=float(y[2]),inward_CN_mg_kg=float(incoming),escaped_CN_mg_kg=float(y[3]),mobile_CN_mg_kg=float(sum(y[:2])),mass_error_mg_kg=float(error))
rows=[]
for n in [1.,1.5,2.]:
 k=math.exp(brentq(lambda logk:integrate(math.exp(logk),n,24.75,1,1)['retained_CN_mg_kg']-ref,-30,0,xtol=1e-11))
 reference=integrate(k,n,24.75,1,1)
 assert abs(reference['retained_CN_mg_kg']-ref)<1e-5
 entry=dict(n=n,k=k,reference=reference,scenarios=[])
 for label,duration,cycles,ironfactor in [('short',.24,400,1),('judgment_midpoint',.5,400,1),('long',1.2,400,1),('R12_delousing',6,270,.85),('R13_delousing',6,270,.9)]:
  r=integrate(k,n,duration,cycles,g)
  r.update(label=label,duration_hours=duration,cycles=cycles,iron_factor=ironfactor,reported_basis_retained_CN_mg_kg=ironfactor*r['retained_CN_mg_kg'])
  entry['scenarios'].append(r)
 rows.append(entry)
output=dict(Q_CN_mg_kg=q,reference_CN_mg_kg=ref,gas_ratio=g,primary_judgment='n=1 as the simplest transfer baseline;0.5h is a rounded intermediate equivalent exposure, not measured contact time. n=1.5 and2 are nonlinear alternatives, not assigned probabilities. No fit to R3,R12,R13.',rows=rows)
Path('educated_estimate_results.json').write_text(json.dumps(output,indent=2)+'\n')
print(json.dumps(output,indent=2))