program xyz2geodetic
    implicit none

    ! Precision double precision (8 Bytes)
    integer, parameter :: dp = kind(1.0d0)
    real(dp), parameter :: PI = 3.14159265358979323846264338327950288_dp
    real(dp), parameter :: RAD2DEG = 180.0_dp / PI

    ! GRS80 Ellipsoid-Parameter (ETRS89)
    real(dp), parameter :: a = 6378137.0_dp                  ! Große Halbachse [m]
    real(dp), parameter :: f = 1.0_dp / 298.257222101_dp     ! Abplattung
    real(dp) :: b, e2, ep2

    ! Eingabekoordinaten (Beispiel aus dem Bild)
    real(dp) :: X, Y, Z

    ! Berechnete geodätische Koordinaten
    real(dp) :: lat, lon, h

    ! Hilfsgrößen
    real(dp) :: p, theta, N

    ! Ellipsoid-Konstanten ableiten
    b = a * (1.0_dp - f)
    e2 = (a**2 - b**2) / (a**2)      ! Erste Exzentrizität im Quadrat
    ep2 = (a**2 - b**2) / (b**2)     ! Zweite Exzentrizität im Quadrat

    ! Testwerte setzen (Geodätischer Referenzpunkt)
    X = 4175739.70_dp
    Y =  853475.32_dp
    Z = 4729951.00_dp

    ! 1. Längengrad lambda
    lon = atan2(Y, X)

    ! 2. Hilfsparameter p und theta (Bowring-Algorithmus)
    p = sqrt(X**2 + Y**2)
    theta = atan2(Z * a, p * b)

    ! 3. Breitengrad phi
    lat = atan2(Z + ep2 * b * (sin(theta)**3), p - e2 * a * (cos(theta)**3))

    ! 4. Ellipsoidische Höhe h
    N = a / sqrt(1.0_dp - e2 * (sin(lat)**2))
    h = (p / cos(lat)) - N

    ! Ausgabe
    print '(A)', "=== ETRS89 XYZ -> Geographische Koordinaten ==="
    print '(A, F12.2, A, F12.2, A, F12.2, A)', "Eingabe:  X = ", X, " m, Y = ", Y, " m, Z = ", Z, " m"
    print '(A)', "------------------------------------------------"
    print '(A, F14.8, A)', "Breite (lat): ", lat * RAD2DEG, " deg N"
    print '(A, F14.8, A)', "Länge  (lon): ", lon * RAD2DEG, " deg E"
    print '(A, F10.2, A)', "Höhe   (h)  : ", h, " m"

end program xyz2geodetic