program ut_correction
  implicit none

  ! Präzision: Doppelte Genauigkeit (64 Bit, C-double)
  integer, parameter :: dp = selected_real_kind(15, 307)

  ! Datumsvariablen
  integer :: year, month, day
  integer :: date_time_vals(8)
  integer :: num_args, stat
  character(len=32) :: arg1, arg2, arg3

  ! Mathematische & Astronomische Variablen
  real(dp) :: pi, jd, mjd, t_bessel
  real(dp) :: delta_ut2_ut1, ut1_utc
  real(dp) :: x_pole, y_pole
  real(dp) :: angle_a, angle_c
  real(dp) :: dt_mjd, s_xy, s_t
  real(dp), parameter :: tai_utc = 37.0_dp

  ! 1. Mathematische Konstante pi
  pi = 4.0_dp * atan(1.0_dp)

  ! 2. Argumente der Kommandozeile auswerten
  num_args = command_argument_count()

  select case (num_args)
  case (0)
    ! Kein Argument -> Aktuelles Systemdatum abfragen
    call date_and_time(values=date_time_vals)
    year  = date_time_vals(1)
    month = date_time_vals(2)
    day   = date_time_vals(3)

  case (1)
    ! Ein Argument (Erwartetes Format: YYYY-MM-DD)
    call get_command_argument(1, arg1)
    if (len_trim(arg1) == 10 .and. arg1(5:5) == '-' .and. arg1(8:8) == '-') then
      read(arg1(1:4), *, iostat=stat) year
      read(arg1(6:7), *, iostat=stat) month
      read(arg1(9:10), *, iostat=stat) day
    else
      print *, "Fehler: Format 'YYYY-MM-DD' erwartet. Beispiel: ./ut_calc 2026-09-06"
      stop
    end if

  case (3)
    ! Drei Argumente (Erwartet: YYYY MM DD)
    call get_command_argument(1, arg1)
    call get_command_argument(2, arg2)
    call get_command_argument(3, arg3)
    
    read(arg1, *, iostat=stat) year
    read(arg2, *, iostat=stat) month
    read(arg3, *, iostat=stat) day

  case default
    print *, "Nutzung:"
    print *, "  Systemdatum verwenden  : ./ut_calc"
    print *, "  Einzelnes Argument     : ./ut_calc YYYY-MM-DD"
    print *, "  Drei Argumente         : ./ut_calc YYYY MM DD"
    stop
  end select

  ! 3. Julianisches Datum (JD) und Modifiziertes Julianisches Datum (MJD)
  jd  = julian_day(year, month, day)
  mjd = jd - 2400000.5_dp

  ! 4. Datum in Bessel-Jahre (T) umrechnen
  t_bessel = 1900.0_dp + (jd - 2415020.31352_dp) / 365.242198781_dp

  ! 5. UT2 - UT1 Korrektur (Saisonale Erddrehungsschwankung)
  delta_ut2_ut1 = 0.022_dp * sin(2.0_dp * pi * t_bessel) - 0.012_dp * cos(2.0_dp * pi * t_bessel) &
                - 0.006_dp * sin(4.0_dp * pi * t_bessel) + 0.007_dp * cos(4.0_dp * pi * t_bessel)

  ! 6. Hilfswinkel A und C für Polbewegung
  angle_a = 2.0_dp * pi * (mjd - 61286.0_dp) / 365.25_dp
  angle_c = 2.0_dp * pi * (mjd - 61286.0_dp) / 435.0_dp

  ! 7. Polkoordinaten x und y [Bogensekunden]
  x_pole =  0.1509_dp + 0.1109_dp * cos(angle_a) - 0.0485_dp * sin(angle_a) &
           - 0.0610_dp * cos(angle_c) - 0.0072_dp * sin(angle_c)

  y_pole =  0.3790_dp - 0.0386_dp * cos(angle_a) - 0.1027_dp * sin(angle_a) &
           - 0.0072_dp * cos(angle_c) + 0.0610_dp * sin(angle_c)

  ! 8. UT1 - UTC Berechnen
  ut1_utc = -0.0737_dp - 0.00024_dp * (mjd - 61294.0_dp) - delta_ut2_ut1

  ! 9. Genauigkeitsabschätzung (S_xy und S_t)
  dt_mjd = mjd - 61286.0_dp
  if (dt_mjd > 0.0_dp) then
    s_xy = 0.00068_dp * (dt_mjd ** 0.80_dp)
    s_t  = 0.00025_dp * (dt_mjd ** 0.75_dp)
  else
    s_xy = 0.0_dp
    s_t  = 0.0_dp
  end if

  ! 10. Ausgabe der Ergebnisse
  print '(A)', "=================================================="
  print '(A)', " ERDROTATIONS- UND POLBEWEGUNGS-PRÄDIKTION"
  print '(A)', "=================================================="
  print '(A, I4.4, A, I2.2, A, I2.2)', " Datum          : ", year, "-", month, "-", day
  print '(A, F14.4)',                  " Julian Day (JD): ", jd
  print '(A, F14.4)',                  " MJD            : ", mjd
  print '(A, F14.6, A)',               " Bessel-Jahr (T): ", t_bessel, " B"
  print '(A)', "--------------------------------------------------"
  print '(A, F10.6, A)',               " UT2 - UT1      : ", delta_ut2_ut1, " s"
  print '(A, F10.6, A)',               " UT1 - UTC      : ", ut1_utc,       " s"
  print '(A, F10.1, A)',               " TAI - UTC      : ", tai_utc,        " s"
  print '(A)', "--------------------------------------------------"
  print '(A, F10.4, A)',               " Polkoordinate x: ", x_pole, " '"
  print '(A, F10.4, A)',               " Polkoordinate y: ", y_pole, " '"
  print '(A)', "--------------------------------------------------"
  print '(A, F10.4, A)',               " Unsicherheit S_xy : ", s_xy, " '"
  print '(A, F10.4, A)',               " Unsicherheit S_t  : ", s_t,  " s"
  print '(A)', "=================================================="

contains

  pure function julian_day(y, m, d) result(jd_val)
    integer, intent(in) :: y, m, d
    real(dp) :: jd_val
    integer :: a, b, c, d_calc

    a = (14 - m) / 12
    b = y + 4800 - a
    c = m + 12 * a - 3
    
    d_calc = d + (153 * c + 2) / 5 + 365 * b + b / 4 - b / 100 + b / 400 - 32045
    
    jd_val = real(d_calc, dp) - 0.5_dp
  end function julian_day

end program ut_correction