from functools import reduce
import argparse
import logging

import cantera as ct

logger = logging.getLogger()
stream_handler = logging.StreamHandler()
logger.addHandler(stream_handler)

parser = argparse.ArgumentParser(description='Calculate Thermodynamic States of the Coal fired Boiler')
parser.add_argument('--hhv', action='store_true', help='Higher heating value is used')
parser.add_argument('-v', '--verbose', action='store_true', help='Verbose output')
args = parser.parse_args()

is_HHV = vars(args)['hhv']
is_verbose = vars(args)['verbose']

if is_verbose:
    logger.setLevel(logging.INFO)
else:
    logger.setLevel(logging.WARNING)

stdT = 298.15       # Temperature at standard state, K

"""#############################################################################
Input Parameter Section
#############################################################################"""

# Air stream Temperatures and Mass flow rates
t = list(map(float, ''' 348.15  315.65  600.79  308.54   318.03  306.11  339.45 '''.split()))
m = list(map(float, ''' 114.4   6.94    362.92  7.25     7.25    7.25    7.25 '''.split()))

coalT = 348.15      # K
coalMfr = 56.8813   # kg/s

# Ultimate Analysis - percentage by weight 
coal_ua = {
    "C":  81.41,
    "H2":  5.47,
    "O2": 10.83,
    "N2": 1.72,
    "S":   0.57,
}

# Proximate Analysis - fraction by weight 
moist   = 0.1551
vm      = 0.3047
fc      = 0.4347
ash     = 0.1055

vm_fc = 1. - moist - ash
fuelMfr = coalMfr

# add mass of moisture and ash
coal_ua['N2'] += 100 * ash / vm_fc
coal_ua["H2O"] = 100 * moist / vm_fc

# Boiler Performance Data
#  - Heat Transfer Rates
heatXferPipe = 677656911.111 # J/s
heatXferWall = 592070866.667 # J/s

# Loading thermo data of species to consider
gri_species = {S.name: S for S in ct.Species.listFromFile('gri30.cti')}
nasa_species = {S.name: S for S in ct.Species.listFromFile('nasa_gas.cti')}
nasa_species.update(gri_species)
system_species = [nasa_species[S] for S in (str.split(' C S H2 O2 N2 CO CO2 H2O SO2 '))]


"""#############################################################################
Calculate averaged air property
#############################################################################"""

gases = [ct.Solution(thermo='IdealGas', species=system_species) for i in range(7)]
for i, gas in enumerate(gases):
    gas.TPY = t[i], ct.one_atm, "O2: 0.233, N2: 0.767"

airs = [
    ct.Quantity(gas, constant='HP', mass=m[i])
    for i, gas in enumerate(gases)
    ]

for i, air in enumerate(airs):
    logger.info("air {}, T = {}, mass flow rate = {}".format(i+1, air.T, air.mass))

airmix = reduce(lambda a, b: a+b, airs)
print("Total Air  flow rate    = ", airmix.mass)


"""#############################################################################
Dummy gaseous coal object containing elementary composition
    mass of ash contents is added to N2
#############################################################################"""

coal = ct.Solution(thermo='IdealGas', species=system_species)
coal.TPY = coalT, ct.one_atm, coal_ua


"""#############################################################################
Coal + v O2 = (Y_C/W_C) CO2 + (Y_H/W_H/2) H2O + (Y_S/W_S) SO2 + (Y_N/W_N) N2
Heating value = Sum(Product Enthalpy of Formation) - Sum(Reactant Enthalpy of Formation)
#############################################################################"""

# coal heating value = 5761 Kcal/kg * 4.184 kJ/kcal
coalHV    = - 5761 * 4.184          # J/kg, with Moisture and Ash
coalHVdaf = - 5761 * 4.184 / vm_fc  # kJ/kg, Dry and Ash Free, Negative since exothermic
print("HHV" if is_HHV else "LHV", "of Coal     , kJ/kg  = ", coalHV)
print("HHV" if is_HHV else "LHV", "of Coal(daf), kJ/kg  = ", coalHVdaf)

# sum product of 1kg coal enthalpy of formation - kJ/kg
def hf_product_coefs (product_name, element_name):
    """ returns kJ/kg """
    # J/kmol / kg/kmol / 1000
    if product_name == "H2O" and is_HHV:
        return (-285820 * 1000.
                / ct.Element(element_name).weight
                / nasa_species[product_name].composition[element_name]
                / 1000. )
    else:
        return (nasa_species[product_name].thermo.h(stdT)
                / ct.Element(element_name).weight
                / nasa_species[product_name].composition[element_name]
                / 1000. )

'''
H(stdT)/W values from NASA polynomial
-  32762.281048240053
- 119952.68929829611
-   9258.666766366705

H_f/W values from google search data
-  32762.45348
- 141887.60121
-   8919.38250

Discrepency in Enthalpy of Formation for H2O is due to phase difference
value above is for vapor and otherwise is for liquid water
'''

logger.info("hf(CO2) / W(C)          = {}".format(hf_product_coefs("CO2", "C")))
logger.info("hf(H2O) / W(H)          = {}".format(hf_product_coefs("H2O", "H")))
logger.info("hf(SO2) / W(S)          = {}".format(hf_product_coefs("SO2", "S")))

sum_product_hf = (
          hf_product_coefs("CO2", "C") * coal.elemental_mass_fraction('C')
        + hf_product_coefs("H2O", "H") * coal.elemental_mass_fraction('H')
        + hf_product_coefs("SO2", "S") * coal.elemental_mass_fraction('S'))

sum_coal_hf = - coalHV + sum_product_hf

logger.info("Sum(Hf_product),  kJ/kg = {}".format(sum_product_hf))
logger.info("Sum(Hf_reactant), kJ/kg = {}".format(sum_coal_hf))


"""#############################################################################
Coal Enthalpy at 348.15 K = \Delta H_f + (H(348.15) - H(298.15))
 - Use graphite's Cp to calculate H difference (H(348.15) - H(298.15))
#############################################################################"""

gr = ct.Solution('graphite.cti')
gr.TP = coalT, ct.one_atm
coal_preheat_enthalpy = gr.enthalpy_mass / 1000.    # kJ/kg
coal_enthalpy = sum_coal_hf + coal_preheat_enthalpy

logger.info("Coal preheat H  , kJ/kg = ", coal_preheat_enthalpy)
logger.info("Coal enthalpy   , kJ/kg = ", coal_enthalpy)
logger.info("Dummy Coal H    , kJ/kg = ", coal.enthalpy_mass/1000.)


"""#############################################################################
Can't Set coal object's H to coal_enthalpy since it is composed of gaseous
    C S H2 O2 N2 and can't have H=coal_enthalpy with T >= 0 K.
Therefore only difference between real coal enthalpy and dummy gas coal is
    caculated and it will be added later.
#############################################################################"""

enthalpy_added_after_mixing = (coal_enthalpy*1000 - coal.enthalpy_mass) * fuelMfr # J
logger.info("enthalpy to add later   = ", enthalpy_added_after_mixing)


################################################################################
## Adiabatic Flame Temperature
##  - also adjust enthalpy difference btw dummy gas fuel and real solid fuel
################################################################################

reactant = airmix + ct.Quantity(coal, constant='HP', mass=fuelMfr)
print("Total Mass flow rate    = ", reactant.mass)


################################################################################
## Adiabatic Flame Temperature
##  - also adjust enthalpy difference btw dummy gas fuel and real solid fuel
################################################################################

reactant.equilibrate("HP")
reactant.HP = reactant.enthalpy/reactant.mass + enthalpy_added_after_mixing/reactant.mass, ct.one_atm
reactant.equilibrate("HP")

print("")
print("Chemical Equilibrium")
reactant.phase()


################################################################################
## Heat Transfer to pipes and walls
################################################################################

reactant.HP = reactant.enthalpy/reactant.mass - (heatXferWall+heatXferPipe)/reactant.mass, ct.one_atm
print("After Heat Transfer")
reactant.phase()

note_ref='''\
Reference for Estimation of Coal Enthalpy of Formation

Eqn (14) and (15) of
Sciazko, M. (2013). Rank-dependent formation enthalpy of coal. Fuel, 114, 2-9. doi:10.1016/j.fuel.2012.06.099
'''