import cantera as ct
import matplotlib.pyplot as plt
import numpy as np

# Inputs from user

spec_sensor = "CH4" # If several species write "CH4-H2"

mixture_name = "Lu_ARC"
is_ARC = True

p = 1e5
tin = 300.0
phi = 1.0
fuel = {'CH4':1}
oxidizer = {'O2': 1, 'N2': 3.76}

# Useful method
def get_delta_T(my_free_flame):
    """
    Compute the thermal thickness of the flame

    author : Q. Douasbin
    modified : T. Lesaffre
    """
    x = my_free_flame.grid
    T = my_free_flame.T
    maxgrad = np.max(np.gradient(T,x))
    delta_T = (max(T)-min(T)) / maxgrad

    return delta_T

# Run flame

## Compile arc mechanism with associated fortran
if is_ARC :
    ct.compile_fortran(f'{mixture_name}.f90')

## Set solution object
gas = ct.Solution(f'{mixture_name}.yaml')
gas.TP = tin, p
gas.set_equivalence_ratio(phi, fuel, oxidizer)

## Set and solve the flame
loglevel = 1 
refine_grid = 'refine'

f = ct.FreeFlame(gas, width=0.02)
f.inlet.X = gas.X
f.inlet.T = gas.T

f.energy_enabled = True
f.set_refine_criteria(ratio=2.0, slope=0.05, curve=0.05, prune=0.01)
f.set_max_jac_age(10, 10)
f.set_time_step(1e-8, [10, 20, 40, 80, 100, 100, 150])
f.max_time_step_count = 500

f.solve(loglevel, refine_grid)

## Get useful information for sensors

sl = f.velocity[0]

thermal_thick = get_delta_T(f)

Omega = np.zeros(f.flame.n_points, 'd')

for spec in spec_sensor.split("-"):
    i_spec = gas.species_index(spec)
    Omega[:] += f.net_production_rates[i_spec, :] \
                        * gas.molecular_weights[i_spec]

print("############################################")
print(f"\nMaximum omega 1D flame based on {spec_sensor}:")
print(float("{:.8g}".format(np.max(abs(Omega)))))
print(f"\nLaminar flame speed:")
print(float("{:.8g}".format(sl)))
print(f"\nThermal flame thickness:")
print(float("{:.8g}".format(thermal_thick)))
print("\n")