  We present fast algorithms to automatically estimate the statistical position
of the bow shock from spacecraft data, using existing analytical
two-dimensional (2D) and three-dimensional (3D) models of the shock surface. We
derive expressions of the standoff distances in 2D and 3D and of the normal to
the bow shock at any given point on it.
  Two simple bow shock detection algorithms are constructed, one solely based
on a geometrical predictor from existing models, the other using this predicted
position to further refine it with the help of magnetometer data, an instrument
flown on many planetary missions.
  Both empirical techniques are applicable to any planetary environment with a
defined shock structure. Applied to the Martian environment and the NASA/MAVEN
mission, the predicted shock position is on average within 0.15 planetary
radius $R_p$ of the bow shock crossing. Using the predictor-corrector
algorithm, this estimate is further refined to within a few minutes of the true
crossing $\approx$0.05 $R_p$).
  Between 2014 and 2021, we detect 14,929 clear bow shock crossings,
predominantly quasi-perpendicular. Thanks to 2D conic and 3D quadratic fits, we
investigate the variability of the shock surface with respect to Mars Years
(MY), solar longitude (Ls) and solar EUV flux levels. Although asymmetry in $Y$
and $Z$ Mars Solar Orbital coordinates is on average small, we show that for
MY32 and MY35, Ls = [135-225$^\circ$] and high solar flux, it can become
particularly noticeable, and is superimposed to the usual North-South asymmetry
due in part to the presence of crustal magnetic fields.
