  In this work we revisit the problem of the dynamical stability of
hierarchical triple systems with applications to circumbinary planetary orbits.
We carry out more than 3 10^8 numerical simulations of planets between the size
of Mercury and the lower fusion boundary (13 Jupiter masses) which revolve
around the center of mass of a stellar binary over long timescales. For the
first time, three dimensional and eccentric planetary orbits are considered. We
explore systems with a variety of binary and planetary mass ratios, binary and
planetary eccentricities from 0 to 0.9 and orbital mutual inclinations ranging
from 0 to 180 degrees. The simulation time is set to 10^6 planetary orbital
periods. We classify the results of those long term numerical integrations into
three categories: stable, unstable and mixed. We provide empirical expressions
in the form of multidimensional, parameterized fits for the two borders that
separate the three dynamical domains . In addition, we train a machine learning
model on our data set in order to have an alternative tool of predicting the
stability of circumbinary planets. Both the empirical fits and the machine
learning model are tested against randomly generated circumbinary systems with
very good results regarding the predictions of orbital stability. The empirical
formulae are also applied to the Kepler and TESS circumbinary systems,
confirming the stability of the planets in these systems. Finally, we present a
REST API with a web based application for convenient access of our simulation
data set.
