How to use the chempy.Reaction function in chempy

def get_rsys(Af=1e16, Ab=1.5e15, Ea=72e3, Er=-12e3):
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, ArrheniusParam(Af, Ea))
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, ArrheniusParam(Ab, Ea-Er))
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

def get_rsys(kf, kb):
    rf = MassAction([kf], unique_keys=['kf'])
    rb = MassAction([kb], unique_keys=['kb'])
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, rf)
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, rb)
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

def get_rsys(Af=1e16, Ab=1.5e15, Ea=72e3, Er=-12e3):
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, ArrheniusParam(Af, Ea))
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, ArrheniusParam(Ab, Ea-Er))
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

def get_rsys(kf, kb):
    rf = MassAction([kf], unique_keys=['kf'])
    rb = MassAction([kb], unique_keys=['kb'])
    fw = Reaction({'Fe+3', 'SCN-'}, {'FeSCN+2'}, rf)
    bw = Reaction({'FeSCN+2'}, {'Fe+3', 'SCN-'}, rb)
    return ReactionSystem([fw, bw], 'Fe+3 SCN- FeSCN+2'.split())

def get_rsys():
    r1 = Reaction({'A'}, {'B'}, MassAction([4./100], unique_keys=['k1']), name='R1: A cons.')
    r2 = Reaction({'B', 'C'}, {'A', 'C'}, MassAction([1e4], unique_keys=['k2']), name='R2: A reform.')
    r3 = Reaction({'B': 2}, {'B', 'C'}, MassAction([3e7], unique_keys=['k3']), name='R3: C form.')
    return ReactionSystem([r1, r2, r3])

def get_rsys():
    r1 = Reaction({'A'}, {'B'}, MassAction([4./100], unique_keys=['k1']), name='R1: A cons.')
    r2 = Reaction({'B', 'C'}, {'A', 'C'}, MassAction([1e4], unique_keys=['k2']), name='R2: A reform.')
    r3 = Reaction({'B': 2}, {'B', 'C'}, MassAction([3e7], unique_keys=['k3']), name='R3: C form.')
    return ReactionSystem([r1, r2, r3])

variables=None):
    pk, = r.param.unique_keys
    if isinstance(r.param, MassAction):
        ratcoeff = to_unitless(
            p[pk], unit_conc**(1-r.order())/unit_time
        )
        if not r.inact_reac:
            r_str = '{}, {}'.format(parmap[pk], r.string(substances=substmap, with_param=False,
                                                         Reaction_arrow='-->', Reaction_coeff_space=''))
        else:
            all_keys = r.keys()
            reac_stoichs = r.all_reac_stoich(all_keys)
            act_stoichs = r.active_reac_stoich(all_keys)
            rate = '*'.join([parmap[pk]] + [('%s^%d' % (substmap[k], v)) if v > 1 else substmap[k]
                                            for k, v in zip(all_keys, act_stoichs) if v > 0])
            r2 = Reaction(dict([(k, v) for k, v in zip(all_keys, reac_stoichs) if v]), r.prod)
            r_str = '{}, {}'.format(rate, r2.string(substances=substmap, with_param=False,
                                                    Reaction_arrow='\u21D2', Reaction_coeff_space=''))
    elif isinstance(r.param, RadiolyticBase):
        ratcoeff = to_unitless(
            p[pk]*variables['doserate']*variables['density'],
            unit_conc/unit_time
        )
        assert not r.reac and not r.inact_reac and not r.inact_prod
        (prod, n), = r.prod.items()
        assert n == 1
        r_str = ('{}, 0 \u21D2 {}' if ratcoeff > 0 else '{}, {} \u21D2 0').format(
            parmap[pk], substmap[prod])
    else:
        raise NotImplementedError("Whats that?")
    return r_str, pk, abs(ratcoeff)

def get_rsys():
    r1 = Reaction({'A'}, {'B'}, MassAction([4./100], unique_keys=['k1']), name='R1: A cons.')
    r2 = Reaction({'B', 'C'}, {'A', 'C'}, MassAction([1e4], unique_keys=['k2']), name='R2: A reform.')
    r3 = Reaction({'B': 2}, {'B', 'C'}, MassAction([3e7], unique_keys=['k3']), name='R3: C form.')
    return ReactionSystem([r1, r2, r3])