  By the virtues of the Dyson-Schwinger equations, we upgrade the published
code \mtt{HELAC} to be capable to calculate the heavy quarkonium helicity
amplitudes in the framework of NRQCD factorization, which we dub
\mtt{HELAC-Onia}. We rewrote the original \mtt{HELAC} to make the new program
be able to calculate helicity amplitudes of multi P-wave quarkonium states
production at hadron colliders and electron-positron colliders by including new
P-wave off-shell currents. Therefore, besides the high efficiencies in
computation of multi-leg processes within the Standard Model, \mtt{HELAC-Onia}
is also sufficiently numerical stable in dealing with P-wave quarkonia (e.g.
$h_{c,b},\chi_{c,b}$) and P-wave color-octet intermediate states. To the best
of our knowledge, it is a first general-purpose automatic quarkonium matrix
elements generator based on recursion relations on the market.
