How to recreate coordinate transformation from vesselfinder (EPSG:4326 to EPSG:3857)?

Question

I'm currently trying to figure out how vesselfinder.com calculates its Box Boundaries (bbox) which they use to query data from their backend.

Given an input like: lat, lon = 59.8230, 22.9586

They fetch data by using this bbox: 13761899,35886447,13779795,35898097

If I try to get a similar bbox by using bboxfinder.com, I get the following values, which aren't even close to what I was expecting: 2553560.4710,8358928.9331,2556565.4293,8360514.8411

The website above is using EPSG:4326 (WGS 84) to EPSG:3857 (WHS 84 / Pseudo-Mercator) by default. I tried to verify in the JS code of vesselfinder that they're using this conversion as well.

    var c = new s.geom.MultiLineString(t);
    return c.transform('EPSG:4326', 'EPSG:3857'),

There are also the following ones mentioned, but I'm pretty sure, that it has to be the upper shown transformation.

it = [
  new $('EPSG:3857'),
  new $('EPSG:102100'),
  new $('EPSG:102113'),
  new $('EPSG:900913'),

The questions now are: What am I doing wrong? / Where do I think wrong?

I also tried using Python for the conversion and even tried out the other mentioned EPSG:XXXXXX types, but haven't got the desired result. I also changed the order of both EPSG types when creating the Transformer, but again, not the desired results.

from pyproj import Transformer

TRAN_4326_TO_3857 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

lat = 59.823002
lon = 22.958583
expansion = 2000

res = TRAN_4326_TO_3857.transform(lng, lat)
bbox = (round(res[0]-expansion), round(res[1]-expansion), round(res[0]+expansion), round(res[1]+expansion))

print(bbox)
# (2455738, 8260436, 2655738, 8460436)

This one is close to the one I got from bboxfinder, but is again not even close to the bbox vesselfinder is using.

Answer 1

https://gis.stackexchange.com/a/370496 seems to have the math.

convertCoordinates(lon, lat) {
    var x = (lon * 20037508.34) / 180;
    var y = Math.log(Math.tan(((90 + lat) * Math.PI) / 360)) /
                  (Math.PI / 180);
    y = (y * 20037508.34) / 180;
    return [x, y];
}

Or, in C# ( https://gis.stackexchange.com/a/325551 )

public static double[] latLonToMeters(double lat, double lon)
{
    //Debug.Log(Mathd.Tan((90d + lat) * Mathd.PI / 720));
    //Debug.Log(Mathd.Tan((90d + lat) * Mathd.PI / 360d));
    double[] meters = new double[2];
    meters[0] = lat * originShift / 180d;
    meters[1] = Mathd.Log(Mathd.Tan((90d+lon) * Mathd.PI / 360d)) /
                   (Mathd.PI / 180d);
    //meters[1] = Mathd.Log(5) / (Mathd.PI / 180d);
    meters[1] = meters[1] * originShift / 180d;
    return meters;
}

In any case, note the website those came from; that may be a better place to get the algorithm. (Then if you need help converting to your preferred language, come back here.)