"""Frozen outputs of the independent Amesim 2404 helium property library."""
import json
import math
from pathlib import Path
import subprocess
import tempfile
import unittest

from app.simulation.native_codegen.build import toolchain

ROOT = Path(__file__).resolve().parents[1]
REFERENCE = ROOT / 'tests/data/amesim_helium_high_temperature_properties.json'


def probe_source(rows, properties=None):
    properties = properties or ROOT / 'native/components/modules/properties.c'
    pairs = ',\n'.join('{%.17g,%.17g,%.17g}' % (r['p'], r['T'], r['h']) for r in rows)
    # The old snapshot is used only by the before/after diagnostic script.
    ideal_cv = 'cp_ideal(T)-rg' if 'cp_ideal(' in properties.read_text() else '1.5*rg'
    return '#include "%s"\n' % properties.as_posix() + r'''
#include <stdio.h>
int main(void) {
    const double inputs[][3]={PAIRS};
    for(unsigned i=0;i<sizeof(inputs)/sizeof(inputs[0]);i++) {
        double p=inputs[i][0],T=inputs[i][1],rho=density(p,T),a,da,dda;
        double v=MOLAR_MASS/rho;attraction(T,&a,&da,&dda);
        double d=v*(v+pr_b)+pr_b*(v-pr_b),dpT=RU/(v-pr_b)-da/d;
        double dpR=(-RU*T/pow(v-pr_b,2)+a*2*(v+pr_b)/(d*d))*(-MOLAR_MASS/(rho*rho));
        double cv=IDEAL_CV+T*dda*log_volume(rho),cp=cv+T*dpT*dpT/(rho*rho*dpR);
        double h=h_ideal(T)+h_departure(p,T),mu=native_viscosity(&helium_medium,T,0);
        double y[2];NativeGas gas;
        if(!native_medium_init(&helium_medium,p,T,.01,0,y) ||
           !native_medium_gas(&helium_medium,y[0],y[1],.01,&gas))return 2;
        double inv=native_temperature_ph(&helium_medium,p,inputs[i][2]);
        double step=.01;
        double dh=(h_ideal(T+step)+h_departure(p,T+step)-h_ideal(T-step)-h_departure(p,T-step))/(2*step);
        double du=(u_ideal(T+step)+u_departure(T+step,rho)-u_ideal(T-step)-u_departure(T-step,rho))/(2*step);
        NativePropertyState cache_states[8],scratch;NativePropertyCache cache;
        native_properties_init(&cache,cache_states,8);
        NativePropertyState *s=property_pt(&cache,&helium_medium,p,T,&scratch);
        local_isentropic(s);
        NativeGas cached;
        if(!native_medium_gas_context(&cache,&helium_medium,y[0],y[1],.01,&cached))return 3;
        double cached_T=native_temperature_ph_context(&cache,&helium_medium,cached.p,cached.h);
        printf("%.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g %.17g\n",
            p,T,rho,h,cp,cv,mu,1/(rho*cv),dpT/(rho*cv),gas.T,gas.p,inv,dh,du,
            gas.h-gas.u-gas.p/gas.rho,s->isentropic_factor,cached_T,gas.h);
    }
    return 0;
}
'''.replace('PAIRS', pairs).replace('IDEAL_CV', ideal_cv)


FIELDS = 'p T rho h cp cv mu dT dP stateT stateP inverseT dhdT dudT enthalpyIdentity isentropicFactor cachedT stateH'.split()


class NativeHeliumHighTemperatureTests(unittest.TestCase):
    @classmethod
    def setUpClass(cls):
        cls.reference = json.loads(REFERENCE.read_text())['rows']
        try:
            compiler = toolchain()[0]
        except (OSError, RuntimeError, subprocess.SubprocessError) as exc:
            raise unittest.SkipTest(f'Native toolchain unavailable: {exc}')
        with tempfile.TemporaryDirectory(prefix='native-helium-high-T-') as tmp:
            directory = Path(tmp); source = directory/'probe.c'; exe = directory/'probe.exe'
            source.write_text(probe_source(cls.reference))
            build = subprocess.run([compiler,'-std=c11','-O3','-Wall','-Wextra','-Werror',
                '-ffp-contract=off','-fno-fast-math','-static-libgcc','-I',str(ROOT/'native/include'),
                str(source),'-lm','-o',str(exe)],capture_output=True,text=True,timeout=60)
            if build.returncode:
                raise AssertionError(build.stderr)
            run = subprocess.run([str(exe)],capture_output=True,text=True,timeout=30)
            if run.returncode:
                raise AssertionError(run.stderr)
            cls.actual = [dict(zip(FIELDS,map(float,line.split()))) for line in run.stdout.splitlines()]
            if len(cls.actual) != len(cls.reference):
                raise AssertionError(f'Expected {len(cls.reference)} rows, got {len(cls.actual)}')

    def test_same_pressure_temperature_against_amesim(self):
        for expected, actual in zip(self.reference,self.actual):
            with self.subTest(p=expected['p'],T=expected['T']):
                self.assertEqual((actual['p'],actual['T']),(expected['p'],expected['T']))
                for field in ('rho','h','cp','cv','mu','dT','dP'):
                    self.assertTrue(math.isclose(actual[field],expected[field],rel_tol=2e-10,abs_tol=1e-12),
                        (field,actual[field],expected[field]))

    def test_state_reconstruction_and_ph_inverse(self):
        for expected, actual in zip(self.reference,self.actual):
            with self.subTest(p=expected['p'],T=expected['T']):
                for field in ('stateT','inverseT','cachedT'):
                    self.assertAlmostEqual(actual[field],expected['T'],delta=2e-6)
                self.assertAlmostEqual(actual['stateP'],expected['p'],delta=expected['p']*2e-10)
                self.assertAlmostEqual(actual['stateH'],expected['h'],delta=max(abs(expected['h'])*2e-10,1e-6))
                self.assertAlmostEqual(actual['enthalpyIdentity'],0,delta=1e-7)

    def test_caloric_derivatives_and_isentropic_consumers(self):
        for expected, actual in zip(self.reference,self.actual):
            with self.subTest(p=expected['p'],T=expected['T']):
                self.assertAlmostEqual(actual['dhdT'],expected['cp'],delta=expected['cp']*2e-7)
                self.assertAlmostEqual(actual['dudT'],expected['cv'],delta=expected['cv']*2e-7)
                # Recover the EOS derivatives independently from rho(p,T)
                # oracle cp/cv and unit heat-input chamber responses.
                dpT=expected['dP']/expected['dT']
                dpR=expected['T']*dpT*dpT/(expected['rho']**2*(expected['cp']-expected['cv']))
                factor=expected['p']/(expected['rho']*dpR*(expected['cp']/expected['cv']))
                self.assertAlmostEqual(actual['isentropicFactor'],factor,delta=abs(factor)*2e-10)


if __name__ == '__main__':
    unittest.main()