!===============================================================================
! Program: etrs89_to_xyz
! Purpose: Umrechnung von geodaetischen Koordinaten (Breite, Laenge, Hoehe)
!          im ETRS89 / GRS80 Ellipsoid in kartesische Koordinaten (X, Y, Z)
!          bezogen auf den Erdschwerpunkt und die Erdrotationsachse (ECEF).
!===============================================================================
program etrs89_to_xyz
    implicit none

    ! Precision definition (double precision / 64-bit real)
    integer, parameter :: dp = selected_real_kind(15, 30)

    ! Constants for GRS80 Ellipsoid
    real(dp), parameter :: a  = 6378137.0_dp          ! Aequatorradius / Grosse Halbachse [m]
    real(dp), parameter :: f  = 1.0_dp / 298.257222101_dp ! Abplattung
    real(dp), parameter :: pi = 3.141592653589793238462643383279502884_dp

    ! Calculated ellipsoid parameters
    real(dp) :: e2  ! Erste numerische Exzentrizitaet zum Quadrat

    ! Input variables
    real(dp) :: lat_deg, lon_deg, h
    real(dp) :: lat_rad, lon_rad

    ! Intermediate and output variables
    real(dp) :: N   ! Querkruemmungshalbmesser (Prime Vertical Radius)
    real(dp) :: X, Y, Z

    ! Calculate square of first numerical eccentricity (e^2 = 2f - f^2)
    e2 = 2.0_dp * f - f**2

    print *, '========================================================='
    print *, '   ETRS89 / GRS80 Geodetic to ECEF (X, Y, Z) Converter   '
    print *, '========================================================='
    print *

    ! Interactive input
    print *, 'Bitte Breitengrad (Latitude) in Grad eingeben (z.B. 48.173333):'
    read(*, *) lat_deg

    print *, 'Bitte Laengengrad (Longitude) in Grad eingeben (z.B. 11.635):'
    read(*, *) lon_deg

    print *, 'Bitte elliptische Hoehe h in Metern eingeben (z.B. 560.0):'
    read(*, *) h

    ! Convert degrees to radians
    lat_rad = lat_deg * pi / 180.0_dp
    lon_rad = lon_deg * pi / 180.0_dp

    ! Compute Prime Vertical Radius of Curvature (Querkruemmungshalbmesser N)
    ! N = a / sqrt(1 - e^2 * sin^2(phi))
    N = a / sqrt(1.0_dp - e2 * (sin(lat_rad)**2))

    ! Compute ECEF Cartesian Coordinates
    ! X = (N + h) * cos(phi) * cos(lambda)
    ! Y = (N + h) * cos(phi) * sin(lambda)
    ! Z = (N * (1 - e^2) + h) * sin(phi)
    X = (N + h) * cos(lat_rad) * cos(lon_rad)
    Y = (N + h) * cos(lat_rad) * sin(lon_rad)
    Z = (N * (1.0_dp - e2) + h) * sin(lat_rad)

    ! Print Results
    print *
    print *, '---------------------------------------------------------'
    print *, 'ERGEBNISSE / RESULTS:'
    print *, '---------------------------------------------------------'
    write(*, '(A, F15.4, A)') ' X = ', X, ' m'
    write(*, '(A, F15.4, A)') ' Y = ', Y, ' m'
    write(*, '(A, F15.4, A)') ' Z = ', Z, ' m'
    print *, '---------------------------------------------------------'

end program etrs89_to_xyz