Hallo,

ich habe ein Problem bei einer Koordinatentransformation.
Ich habe Koordinaten im GaussKrüger Format (Rechtswert/Hochwert) im Quellkoordinatensystem Pulkovo 1942 (Krassovsky). Ich würde diese gerne ins WGS84 transformieren.

Zum Vorgehen:
1.) Umwandeln von GaussKrüger Koordinaten in Geo Koordinaten
2.) Umwandeln von Geographischen Koordinaten in Geozentrische Koordinaten
3.) Speziallfall Molodensky-Badekas Tranformation
4.) Umwandeln von Geozentrischen Koordinaten in Geographische Koordinaten

Ich nehme als Referenz folgende Online-Transformation:
https://epsg.io/transform#s_srs=2399-1675&t_srs=4326&x=5389306.9210650&y=5823166.2710200

Meine Daten:
Rechtswert: 5389306.921065
Hochwert: 5823166.27102

Was ich erhalten möchte:
Längengrad: 13.3672138°
Breitengrad: 52.5249447°

1.)

        a = 6378245;
        f = 1/293.3;
        c = a / (1 - f);

        // Quadrat der zweiten numerischen Exzentrizität
        ex2 = (2 * f - f * f) / ((1 - f) * (1 - f));
        ex4 = ex2 * ex2;
        ex6 = ex4 * ex2;
        ex8 = ex4 * ex4;

        // Koeffizienten zur Berechnung der geographischen Breite aus gegebener
        // Meridianbogenlänge
        e0 = c * (pi / 180) * (1 - 3 * ex2 / 4 + 45 * ex4 / 64 - 175 * ex6 / 256 + 11025 * ex8 / 16384);
        f2 =            (180 / pi) * (3 * ex2 / 8 - 3 * ex4 / 16 + 213 * ex6 / 2048 - 255 * ex8 / 4096);
        f4 =                          (180 / pi) * (21 * ex4 / 256 - 21 * ex6 / 256 + 533 * ex8 / 8192);
        f6 =                                        (180 / pi) * (151 * ex6 / 6144 - 453 * ex8 / 12288);

        // Geographische Breite bf zur Meridianbogenlänge gf = hw
        sigma = (northing) / e0;
        sigmr = sigma * pi / 180;
        bf = sigma + f2 * Math.sin(2 * sigmr) + f4 * Math.sin(4 * sigmr) + f6 * Math.sin(6 * sigmr);

        // Breite bf in Radianten
        br = bf * pi / 180;
        tan1 = Math.tan(br);
        tan2 = tan1 * tan1;
        tan4 = tan2 * tan2;

        cos1 = Math.cos(br);
        cos2 = cos1 * cos1;

        etasq = ex2 * cos2;

        // Querkrümmungshalbmesser nd
        nd = c / Math.sqrt(1 + etasq);
        nd2 = nd * nd;
        nd4 = nd2 * nd2;
        nd6 = nd4 * nd2;
        nd3 = nd2 * nd;
        nd5 = nd4 * nd;

        // Längendifferenz dl zum Bezugsmeridian lh
        kz = (int)(easting / 1e6);
        lh = kz * 3;
        dy = (easting - (kz * 1e6 + 500000));
        dy2 = dy * dy;
        dy4 = dy2 * dy2;
        dy3 = dy2 * dy;
        dy5 = dy4 * dy;
        dy6 = dy3 * dy3;

        b2 = -tan1 * (1 + etasq) / (2 * nd2);
        b4 = tan1 * (5 + 3 * tan2 + 6 * etasq * (1 - tan2)) / (24 * nd4);
        b6 = -tan1 * (61 + 90 * tan2 + 45 * tan4) / (720 * nd6);

        l1 = 1 / (nd * cos1);
        l3 = -(1 + 2 * tan2 + etasq) / (6 * nd3 * cos1);
        l5 = (5 + 28 * tan2 + 24 * tan4) / (120 * nd5 * cos1);

        // Geographischer Breite bp und Länge lp als Funktion von Rechts- und
        // Hochwert
        latitude = bf + (180 / pi) * (b2 * dy2 + b4 * dy4 + b6 * dy6);
        longitude = lh + (180 / pi) * (l1 * dy + l3 * dy3 + l5 * dy5)


Nach 1.) erhalte ich folgende Werte:

Längengrad: 13.369036382735441
Breitengrad: 52.525467788632824

2.) Umrechnung in geozentrische Koordinaten:

        pi = Math.PI;
        a = 6378245;
        f = 1/293.3;
        e2 = (2 * f - f * f);

        longitudeRAD = longitude * (pi / 180);
        latitudeRAD = latitude * (pi / 180);

        // Querkrümmungshalbmesser
        double nd = a /sqrt(1 - e2 * sin(latitudeRAD)²);
       
        // Umrechnung in geozentrische Koordinaten
        xp = (nd) * cos(latitudeRAD) * cos(longitudeRAD);
        yp = (nd) * cos(latitudeRAD) * sin(longitudeRAD);
        zp = ((1 - e2wgs84) * nd) * sin(latitudeRAD);


Nach 2.) erhalte ich folgende geozentrische Koordinaten:

xp: 3783403.403351258
yp: 899173.4481951335
zp: 5038677.955342111

3.) Transformation von Pulkovo 1942 zu WGS84


Transformationsparameter:
        dx = 24.0;
        dy = -123.0;
        dz = -94.0;
        rx = -0.02 * (pi / 180);
        ry = 0.25 * (pi / 180);
        rz = 0.13 * (pi / 180);
        scale = 1.1;

        M = 1 + scale*10e-6;

        // 7-Parameter Transformation
        x = xp + (M * (dx - (rz * dy) + (ry * dz)));
        y = yp + (M * ((rz * dx) + dy - (rx * dz)));
        z = zp + (M * ((-ry * dx) + (rx * dy) + dz));


4.) Umrechnung geozentrische Koordinaten in geographische Koordinaten


        rb = Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2));
       
        latitutdeWGS84 = (180 / pi) * Math.atan((z / rb) / (1 - e2wgs84));
        longitudeWGS84 = (180 / pi) * Math.atan(y / x);

Hier meine Frage:
Laut https://epsg.io/1032-method multipliziere ich die Quellkoordinaten mit den Rotationen und M und verschiebe dann um dx,dy,dz.
Mache ich es so, erhalte ich die folgenden Werte:

52.27755746434763°, 13.447067545543165°

Wenn ich wie oben in 3.) beschrieben dx,dy und dz rotiere und mit M multiplizieren und dann auf xp,yp,zp addiere erhalte ich folgende Werte:

52.52501509125434°, 13.367192233116334°

Ersteres macht mehr Sinn und ist ja auch so definiert. Die Werte sehen umgedreht allerdings besser aus.

Die unteren Werte haben eine Distanz von 8m zum erwarteten Punkt. Ich hätte es gerne noch genauer und vermute ich habe irgendwo einen Denkfehler.

Ich wäre über jegliche Anregung dankbar!