We present general empirical analytical equations of bow shock structures historically used at Mars, and show how to estimate automatically the statistical position of the bow shock with respect to spacecraft data from 2D polar and 3D quadratic fits. Analytical expressions of bow shock normal in 2D and 3D are given for any point on the shock’s surface. This empirical technique is applicable to any planetary environment with a defined shock structure. Applied to the Martian environment and the NASA/MAVEN mission, the predicted bow shock location from ephemerides data is on average within 0.15 planetary radius R_p of the actual shock crossing as seen from magnetometer data. Using a simple predictor-corrector algorithm based on the absolute median deviation of the total magnetic field and the general form of quasi-perpendicular shock structures, this estimate is further refined to within a few minutes of the true crossing (≈0.05 R_p). With the refined algorithm, 14,929 bow shock crossings, predominantly quasi-perpendicular, are detected between 2014 and 2021. Analytical 2D conic and 3D quadratic surface fits, as well as standoff distances, are successively given for Martian years 32 to 35, for several (seasonal) solar longitude ranges and for two solar EUV flux levels. Although asymmetry in Y and Z Mars Solar Orbital coordinates is on average small, it is shown that for Mars years 32 and 35, Ls = [135-225] degrees and high solar flux, it can become particularly noticeable and is superimposed to the usual North-South asymmetry due to the presence of crustal magnetic fields.