Hallo,

ich stelle Dir einfach mal meinen Matlab-Code zur Verfügung - in der Hoffnung, dass Dein Nickname sich auf die Hochschule Dresden bezieht und Du diesen Code lesen und ggf. ausführen kannst. Zwischenergebnisse folgen weiter unten. Hier erst einmal der Code:

%% Alle (Datums-)Punkte im globalen System
pointsXYZ = [
    3890198.4748    860545.2218     4964144.7203 % P0
    3875828.2460    846780.7380     4977583.1370 % P1
];
 
%% Globaler Fundamentalpunkt
% fundamentalPoint = mean(pointsXYZ);
fundamentalPoint = pointsXYZ(1,:);
X0 = fundamentalPoint(1);          
Y0 = fundamentalPoint(2);
Z0 = fundamentalPoint(3);
 
% Zugehoerige geographische Koordinaten im GRS80
[B0, L0, h0] = XYZ2BLh(X0, Y0, Z0, 'GRS80');
 
%% Lokaler Fundamentalpunkt
x0 = 100000.0;
y0 = 100000.0;
z0 = 100000.0;
 
% Rotationsmatrix, um globale XYZ Koordinaten in lokale Koordinaten zu
% konvertieren, wobei die lokalen Koordinaten in einem geodaetischen ENH
% System definiert sind
 
R = [
             -sind(L0)            cosd(L0)     0
    -sind(B0)*cosd(L0)  -sind(B0)*sind(L0)   cosd(B0)
     cosd(B0)*cosd(L0)   cosd(B0)*sind(L0)   sind(B0)
];
 
 
pointsENH = [x0 y0 z0] + (R * (pointsXYZ - fundamentalPoint)')';
 
fprintf('%32.16f  %32.16f  %32.16f\n', pointsENH');

2. Bestimmung von Länge, Breite und Höhe von P0

λP0 [°] 51,4370585340
ϕP0 [°] 12,4734506530
hP0 [m] 184,8665571660

Ich ganz leicht variierende Werte, wobei dies sicher auf leicht differierende Parameter bei der Definition des Ellipsoidis zurückzuführen ist. Ich nutze a = 6378137.0 und f = 298.257222100880030 für die Berechnung und erhalten dann entsprechend

ϕ0 =  51.4370585336273°
λ0 =  12.4734506531850°
h0 = 184.8665571734310 m

Insofern scheint dies zu passen, bis auf die Bezeichnung λ vs. ϕ. Sofern P0 nicht (näherungsweise) in der Mitte Deines Koordinatensystems liegt, würde ich diesen nicht verwenden. In P0 sind ζx und ζy beide null und werden dann mit zunehmenden Abstand größer.

3. Definition von lokalen Koordinaten von P0, um negative Koordinaten zu vermeiden (optional)

P0: (y)100000,0000 (x)100000,0000 (z)100000,0000

Ich persönlich nutze immer verschiedene Werte für den Ost- und Nordwert, um diese an der Koordinate schon zu erkennen. Die Höhe würde ich nur so riesig wählen, wenn es dafür einen Grund gibt. Beachte bitte, dass auch der Double grenzen hat und Du mit vielen (unnötigen) Vorkommastellen bereits einen (Groß-)Teil der 64 bit belegst. Dies aber mehr als Hinweis.

4. Transformation Deiner geozentrischen Datumskoordinaten in das lokale System mittels P0 und unter Verwendung der Gleichung über Abb. 2.

P0: (Y)100000,0000 (X)100000,0000 (Z)100000,0000
P1: (Y)102656,0032 (X)117380,6603 (Z)83647,2074

Ich erhalte für als Rechts- und Hochwert bzw. Höhe folgende ungerundete Werte für P1

E =  89664.1958921295881737 m
N = 121672.9835989736020565 m
H =  99907.9738623794983141 m

Hier weichen unsere Lösungen ab und ich vermute, Du hast die Drehmatrix nicht korrekt definiert, sodass diese die geozentrischen Werte in der Folge X-Y-Z erwartet, d.h.

$\begin{pmatrix}E_i \\ N_i \\ H_i\end{pmatrix} = \begin{pmatrix} E_0 \\ N_0 \\ H_0 \end{pmatrix} + \begin{pmatrix} -\sin\lambda_0 & \cos\lambda_0 & 0 \\ -\sin\phi_0\cos\lambda_0 & -\sin\phi_0\sin\lambda_0 & \cos\phi_0 \\ \cos\phi_0\cos\lambda_0 & \cos\phi_0\sin\lambda_0 & \sin\phi_0 \end{pmatrix} \begin{pmatrix} X_i - X_0 \\ Y_i - Y_0 \\ Z_i - Z_0 \end{pmatrix}$

Schau mal, ob Du damit weiterkommst Du die Abweichungen findest.

Viele Grüße
Micha

Anhang: Der Code vom Script XYZ2BLh.m lautet:

%% Kartesische X,Y,Z -> Ellipsoidische B,L,h
%   
% @param <Double>X ..... X-Komponente
% @param <Double>Y ..... Y-Komponente
% @param <Double>Z ..... Z-Komponente
% @param <String>ell ... Name des Referenzellipsoid
%
% @return <Double>B .... Breite
% @return <Double>L .... Länge
% @return <Double>h .... Höhe
%
% @version 1.0 vom 30.10.2007 [XYZ2BLh.m]
% @author Michael Loesler
 
function [B, L, h] = XYZ2BLh(X, Y, Z, ell)
    if nargin==3
        ell = 'WGS84';
    elseif nargin<3
        error('Zu wenig Eingangsargumente!');
    end
 
    if strcmp(upper(ell),'WGS84') || strcmp(upper(ell),'GRS80')
        a = 6378137.0;
        alpha = 1.0/298.257222100880030;
    elseif strcmp(upper(ell), 'BESSEL')
        a = 6377397.155;
        alpha = 1.0/299.15281285;
    elseif strcmp(upper(ell), 'KRASSOWSKI')
        a = 6378245.0;
        alpha = 1.0/298.3;
    elseif strcmp(upper(ell), 'HAYFORD')
        a = 6378388.0;
        alpha = 1.0/297.0;
    else
        % Kugel
        a = 6370000;
        alpha = 0;
    end
 
    roh = 180/pi;
 
    b = a*(1-alpha);
    c = a^2/b;
 
    e1 = sqrt((a^2-b^2)/a^2);
    e2 = sqrt((a^2-b^2)/b^2);
 
    p = sqrt(X^2+Y^2);    
 
    Theta = atan((Z*a)/(p*b));
 
    L = atan(Y/X)*roh;
    B = atan((Z+ e2^2*b*(sin(Theta))^3) / (p-e1^2*a*(cos(Theta))^3));
 
    eta2 = e2^2*(cos(B))^2;
    V = sqrt(1 + eta2);
    N = c/V;
 
    h = p/cos(B) - N;
    B = B*roh;
return;