import sys
from math import sqrt
from collections import defaultdict

GAMMA, VOL = 2.3, 37   # C-gamma model parameters

# atomic weights of elements (g/mol)
ATOMW = {'C': 12.011, 'H': 1.008, 'N': 14.007, 'O': 15.999, 'F': 18.998, 'Cl': 35.453,
    'S': 32.067, 'Br': 79.904, 'Na': 22.99, 'Ag': 107.868, 'Al': 26.982, 'Ba': 137.328}

# Formation enthalpies of decomposition products (kJ/mol) from NIST Webbook
FORMATION_ENTHALPY = dict(H2O=-241.826,CO2=-393.51,CO=-110.53, 
    HF=-273.30, HCl=-92.31, CF4=-933.,CCl4=-95.8, HBr=-36.29, 
    Cs=0, N2=0, H2=0, O2=0, F2=0, Cl2=0, Br2=0, S8=0, 
    NaCl=-411.12, Si=0, K=0, PO3=-450, Ag=0, Na=0,
    Al2O3=-1675.7, BaO=-548.10, BaOH=-226.44, BaH2O2=-626.56, Al=0)

GASES = ['CO','CO2','H2O','N2','H2','O2','F2','Cl2','HF','HCl','HBr','Br2']

def str2tup(chunk):  # => used only by str2mf !
    if len(chunk) == 1: return (chunk, 1)
    slen = 2 if chunk[1].isalpha() else 1
    symbol, howmany = chunk[:slen], chunk[slen:]
    if howmany == '': return (symbol, 1)
    return (symbol, eval(howmany.strip()))

def str2mf(string):
    """ read a string as input ans return corresponding formula as dict:
    str2mf('C3H7NO') ==> {'C': 3, 'H': 7, 'N': 1, 'O': 1}
    """
    js = [0]+[j for j in range(1, len(string)) if string[j].isupper()]
    chunks = [string[i:j] for i,j in zip(js, js[1:] + [None])]
    return dict(str2tup(chunk) for chunk in chunks)

def get_prodets(empirical_formula):    # prodHoF = 0
    """ read a dict as input and returns a dict of prodets:
    >>> get_prodets({'C':6,'H':6,'N':6,'O':6})
        ==> {'CO2': 1.5, 'H2O': 3.0, 'Cs': 4.5, 'N2': 3.0}
    """
    if set(empirical_formula) - set(ATOMW): return None  # not implemented elements
    g, p = defaultdict(int, empirical_formula), defaultdict(lambda: 0)
    # Al
    p['Al2O3'] = 0.5*min(g['Al'],(2/3)*g['O'])
    g['O'] -= 3*p['Al2O3']
    g['Al'] -= 2*p['Al2O3']
    # Ba
    p['BaO'] = min(g['O'],g['Ba'])
    g['Ba'] -= p['BaO']
    g['O'] -= p['BaO']
    # Na -> NaCl
    p['NaCl'] = min(g['Na'],g['Cl'])
    g['Na'] -= p['NaCl']
    g['Cl'] -= p['NaCl']
    # HF
    p['HF'] = min(g['H'],g['F'])
    g['H'] -= p['HF']
    g['F'] -= p['HF']
    # CF4
    p['CF4'] = min(g['C'],0.25*g['F'])
    g['C'] -= p['CF4']
    g['F'] -= 4*p['CF4']
    # HCl
    p['HCl'] = min(g['H'],g['Cl'])
    g['H'] -= p['HCl']
    g['Cl'] -= p['HCl']
    # CCl4
    p['CCl4'] = min(g['C'],0.25*g['Cl'])
    g['C'] -= p['CCl4']
    g['Cl'] -= 4*p['CCl4']
    # H2O
    p['H2O'] = min(g['O'],0.5*g['H'])
    g['O'] -= p['H2O']
    g['H'] -= 2*p['H2O']
    # CO2
    p['CO2'] = min(g['C'],0.5*g['O'])
    g['O'] -= 2*p['CO2']
    g['C'] -= p['CO2']
    # CO
    p['CO'] = min(g['O'],g['C'])
    g['C'] -= p['CO']
    g['O'] -= p['CO']
    # HBr
    p['HBr'] = min(g['H'],g['Br'])
    g['H'] -= p['HBr']
    g['Br'] -= p['HBr']
    # single-element products
    p['Cs'] = g['C']
    p['Ag'] = g['Ag']
    p['Na'] = g['Na']
    p['S8'] = 1./8*g['S']
    p['H2'] = 0.5*g['H']
    p['N2'] = 0.5*g['N']
    p['O2'] = 0.5*g['O']
    p['F2'] = 0.5*g['F']
    p['Cl2'] = 0.5*g['Cl']
    p['Br2'] = 0.5*g['Br']
    return p # dict of prodets

def calculate_uG(formula, HoF, TMD, density):
    """ input:
    formula = empirical formula (dict)
    HoF = heat of formation (kJ/mol)
    TMD = Theoretical Maximal Density
    density = loading density
        output:
    Gurney velocity in km/s
    """
    prodets = get_prodets(formula)
    prodHoF = sum(prodets[k]*FORMATION_ENTHALPY[k] for k in prodets)
    gasmol = sum(prodets.get(gas, 0) for gas in GASES)
    molw = sum(o*ATOMW[k] for k,o in formula.items())
    N = gasmol / molw
    Ed = HoF - prodHoF
    Q = Ed / molw
    rho0 = sqrt(density*TMD)
    Vo_by_V = 1/(VOL*N*rho0)
    coef = (GAMMA/(GAMMA+1))**GAMMA
    jem = 1 - 2*coef*Vo_by_V**(GAMMA-1)
    EG = max([jem * Q, 0]) # Q <= 0 implies that the material is NOT energetic
    return sqrt(2*EG)      # Gurney velocity uG


if __name__ == '__main__':
    while True:
        formula = input('Empirical formula = ')
        formula = str2mf(formula)
        try:
            HoF = float(input('Heat of formation (kJ/mol) = '))
            TMD = float(input('Theoretical maximal density (TMD) = '))
            density = input('Loading density (if < TMD) = ')
            density = TMD if density == '' else float(density)
        except ValueError:
            break
        # ready: formula, HoF, TMD, density
        uG = calculate_uG(formula, HoF, TMD, density)
        print(f'>> {uG:5.3f} km/s')