! -------------------------------------------------------------------------
! Bernese GPS Software Version 5.1
! -------------------------------------------------------------------------

MODULE s_ITRF2ETRF

  IMPLICIT NONE

#include "COM_CPU.inc"

CONTAINS

SUBROUTINE itrf2etrf(timcrd,datum,trafo,xstat,xstnew,xvel,xvenew)

! -------------------------------------------------------------------------
! Purpose:    Transformation from ITRFxx to ETRFxx or from ETRFxx to ITRFxx
!
! Author:     S. Lutz
!
! Created:    04-Nov-2011
!
! Changes:    02-Dec-2016 YJ: Add call cpu_enter and cpu_exit
!             10-Apr-2017 ST: Move IMPLICIT NONE to module scope
!             10-Apr-2017 ST: Add #include "COM_CPU.inc"
!             10-Apr-2017 ST: Add IF(srfDebug) in front of CALL cpu_
!             27-MAR-2020 LH: Added [E|I]TRF2014
!             22-Dec-2020 RD: Transform D_CONST to a protected structure CONST
!             28-Nov-2023 SL: Print transformation parameters
!             28-Nov-2023 SL: Support of velocity transformation
!             12-Jan-2024 SL: ITRF2020 to ETRF2020 parameters added
!             31-Jan-2024 SL/SS: Adopting the ETRS89 tool extension made at PNAC
!             01-Feb-2024 SL/SS: Refined calculation of dt
!             12-Mar-2024 SL/SS: Adjusted specification of dt
!
! Copyright:  Astronomical Institute
!             University of Bern
!             Switzerland
! -------------------------------------------------------------------------

! Modules
! -------
  USE m_bern,   ONLY: i4b, r8b, lfnErr, lfnPrt
  USE m_cpu,    ONLY: cpu_enter, cpu_exit
  USE d_const,  ONLY: CONST
  USE f_djul
  USE s_exitrc


! List of arguments
! -----------------
! input:
  REAL(r8b)                          :: timcrd     ! Epoch of coordinates [MJD]
  CHARACTER(LEN=*)                   :: datum      ! Datum string
  INTEGER(i4b)                       :: trafo      ! Sign of tranformation
  REAL(r8b), DIMENSION(3)            :: xstat      ! Input coordinates
  REAL(r8b), DIMENSION(3), OPTIONAL  :: xvel       ! Input velocities

! output:
  REAL(r8b), DIMENSION(3)            :: xstnew     ! Output coordinates
  REAL(r8b), DIMENSION(3), OPTIONAL  :: xvenew     ! Output velocities

! Local Parameters
! ----------------
  CHARACTER(LEN=*), PARAMETER        :: srName = 'itrf2etrf'
  INTEGER(i4b),     PARAMETER        :: ntrf = 12  ! Nbr of supported ITRFs

! Local Variables
! ---------------
  CHARACTER(LEN=2), DIMENSION(ntrf)  :: trflst
  CHARACTER(LEN=2)                   :: yy

  REAL(r8b),        DIMENSION(ntrf,8):: param
  REAL(r8b),        DIMENSION(ntrf,7):: rates
  REAL(r8b),        DIMENSION(8)     :: pari
  REAL(r8b),        DIMENSION(7)     :: rati
  REAL(r8b),        DIMENSION(7)     :: parc, parv
  REAL(r8b),        DIMENSION(3,3)   :: rotc, rotv
  REAL(r8b)                          :: doy, dti

  INTEGER(i4b)                       :: i, j

  LOGICAL,                      SAVE :: first = .TRUE.


! DEBUG: entry
! ------------
  IF (srfDebug) CALL cpu_enter(srName)


! Default values
! --------------
  param(:,:) = 0d0
  rates(:,:) = 0d0

  trflst = (/'89','90','91','92','93','94','96','97','00','05','14','20'/)


! Transformation parameters from ITRFyy to ETRFyy
! -----------------------------------------------
!                  Tx[mm]   Ty[mm]   Tz[mm]    D[ppb]   Rx[mas]   Ry[mas]   Rz[mas]   t0[y]
! ETRF2020
  param(12,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates(12,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.086d0,  0.519d0, -0.753d0/)
! ETRF2014
  param(11,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates(11,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.085d0,  0.531d0, -0.770d0/)
! ETRF2005
  param(10,:) = (/ 56.0d0,  48.0d0, -37.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates(10,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.054d0,  0.518d0, -0.781d0/)
! ETRF2000
  param( 9,:) = (/ 54.0d0,  51.0d0, -48.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 9,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.081d0,  0.490d0, -0.792d0/)
! ETRF97
  param( 8,:) = (/ 41.0d0,  41.0d0, -49.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 8,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.200d0,  0.500d0, -0.650d0/)
! ETRF96
  param( 7,:) = (/ 41.0d0,  41.0d0, -49.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 7,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.200d0,  0.500d0, -0.650d0/)
! ETRF94
  param( 6,:) = (/ 41.0d0,  41.0d0, -49.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 6,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.200d0,  0.500d0, -0.650d0/)
! ETRF93
  param( 5,:) = (/ 19.0d0,  53.0d0, -21.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 5,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.320d0,  0.780d0, -0.670d0/)
! ETRF92
  param( 4,:) = (/ 38.0d0,  40.0d0, -37.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 4,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.210d0,  0.520d0, -0.680d0/)
! ETRF91
  param( 3,:) = (/ 21.0d0,  25.0d0, -37.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 3,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.210d0,  0.520d0, -0.680d0/)
! ETRF90
  param( 2,:) = (/ 19.0d0,  28.0d0, -23.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 2,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.110d0,  0.570d0, -0.710d0/)
! ETRF89
  param( 1,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.000d0,  0.000d0,  0.000d0, 1989.0d0/)
  rates( 1,:) = (/  0.0d0,   0.0d0,   0.0d0,   0.00d0,   0.110d0,  0.570d0, -0.710d0/)


! Get transformation parameters
! -----------------------------
  dti = 0d0
  i = 0

  DO j=1,ntrf
    IF(datum(len_trim(datum)-3:len_trim(datum)) == 'R3_0') THEN
      IF(INDEX(datum(len_trim(datum)-5:len_trim(datum)-4),trflst(j))==0) CYCLE
    ELSEIF(datum(len_trim(datum)-1:len_trim(datum)) == 'R3') THEN
      IF(INDEX(datum(len_trim(datum)-3:len_trim(datum)-2),trflst(j))==0) CYCLE
    ELSEIF(datum(len_trim(datum):len_trim(datum)) == 'P') THEN
      IF(INDEX(datum(len_trim(datum)-2:len_trim(datum)-1),trflst(j))==0) CYCLE
    ELSEIF(datum(len_trim(datum)-1:len_trim(datum)) == '_0') THEN
      IF(INDEX(datum(len_trim(datum)-3:len_trim(datum)-2),trflst(j))==0) CYCLE
    ELSE
      IF(INDEX(datum(len_trim(datum)-2:len_trim(datum)),trflst(j))==0) CYCLE
    ENDIF
    yy   = trflst(j)
    pari = param(j,:)
    rati = rates(j,:)
    doy  = (param(j,8)-REAL(INT(param(j,8))))*365.25d0+1d0
    dti  = (timcrd-DJUL(INT(param(j,8)),1,doy))/365.25d0
    i = i+1
  ENDDO

  IF(i.EQ.0) THEN
    WRITE(lfnErr,"(/,' *** SR ',A,':',A)") srName, &
         ' No transformation parameters found'
    WRITE(lfnErr,'(8X,A,A)') 'Datum: ',trim(datum)
    WRITE(lfnErr,*)
    CALL exitRc(2)
  ENDIF


! Apply sign of transformation
! ----------------------------
  pari(:) = trafo*pari(:)
  rati(:) = trafo*rati(:)

  IF(first) THEN
    WRITE(lfnPrt,*)
    IF(trafo>0) THEN
      WRITE(lfnPrt,"(1X,A,A,A,A,A)") &
           'ITRF-',yy,' to ETRF-',yy,':'
    ELSE
      WRITE(lfnPrt,"(1X,A,A,A,A,A)") &
           'ETRF-',yy,' to ITRF-',yy,':'
    ENDIF
    WRITE(lfnPrt,'(1X,18("-"))')
    WRITE(lfnPrt,*)

    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3A8,A9,3A10,2X,A8)') &
         'Tx','Ty','Tz', &
         'D', &
         'Rx','Ry','Rz', &
         'Epoch'
    WRITE(lfnPrt,'(1X,3A8,A9,3A10,2X,A8)') &
         'mm','mm','mm', &
         'ppb', &
         'mas','mas','mas', &
         'year'
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3F8.2,F9.3,3F10.4,2X,F8.3)') &
         pari(1),pari(2),pari(3), &
         pari(4), &
         pari(5),pari(6),pari(7), &
         ABS(pari(8))
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3A8,A9,3A10)') &
         'mm/y','mm/y','mm/y', &
         'ppb/y', &
         'mas/y','mas/y','mas/y'
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3F8.2,F9.3,3F10.4)') &
         rati(1),rati(2),rati(3), &
         rati(4), &
         rati(5),rati(6),rati(7)
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,*)
  ENDIF

  DO i=1,7
    IF(i==1.OR.i==2.OR.i==3) THEN
      parc(i) = (pari(i)+rati(i)*dti)*1d-3
      parv(i) = rati(i)*1d-3
    ELSEIF(i==4) THEN
      parc(i) = (pari(i)+rati(i)*dti)*1d-9
      parv(i) = rati(i)*1d-9
    ELSEIF(i==5.OR.i==6.OR.i==7) THEN
      parc(i) = (pari(i)+rati(i)*dti)*1d-3/CONST%ars
      parv(i) = rati(i)*1d-3/CONST%ars
    ENDIF
  ENDDO

  IF(first) THEN
    WRITE(lfnPrt,'(1X,A)') &
         'Transformation parameters at epoch of coordinates:'
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3A8,A9,3A10,2X,A8)') &
         'Tx','Ty','Tz', &
         'D', &
         'Rx','Ry','Rz', &
         'DeltaT'
    WRITE(lfnPrt,'(1X,3A8,A9,3A10,2X,A8)') &
         'mm','mm','mm', &
         'ppb', &
         'mas','mas','mas', &
         'years'
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3F8.2,F9.3,3F10.4,2X,F8.3)') &
         parc(1)*1d3,parc(2)*1d3,parc(3)*1d3, &
         parc(4)*1d9, &
         parc(5)*1d3*CONST%ars,parc(6)*1d3*CONST%ars,parc(7)*1d3*CONST%ars, &
         dti
    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,*)

    first = .FALSE.

  ENDIF


! Compute rotation matrix elements
! --------------------------------
  rotc(1,1) =  parc(4)*xstat(1)
  rotc(1,2) = -parc(7)*xstat(2)
  rotc(1,3) =  parc(6)*xstat(3)
  rotc(2,1) =  parc(7)*xstat(1)
  rotc(2,2) =  parc(4)*xstat(2)
  rotc(2,3) = -parc(5)*xstat(3)
  rotc(3,1) = -parc(6)*xstat(1)
  rotc(3,2) =  parc(5)*xstat(2)
  rotc(3,3) =  parc(4)*xstat(3)

  rotv(1,1) =  parv(4)*xstat(1)
  rotv(1,2) = -parv(7)*xstat(2)
  rotv(1,3) =  parv(6)*xstat(3)
  rotv(2,1) =  parv(7)*xstat(1)
  rotv(2,2) =  parv(4)*xstat(2)
  rotv(2,3) = -parv(5)*xstat(3)
  rotv(3,1) = -parv(6)*xstat(1)
  rotv(3,2) =  parv(5)*xstat(2)
  rotv(3,3) =  parv(4)*xstat(3)


! Compute transformation
! ----------------------
  DO i=1,3
    xstnew(i) = xstat(i)
    DO j=1,3
      xstnew(i) = xstnew(i)+rotc(i,j)
    ENDDO
    xstnew(i) = xstnew(i)+parc(i)

    IF(PRESENT(xvel)) THEN
      xvenew(i) = xvel(i)
      DO j=1,3
        xvenew(i) = xvenew(i)+rotv(i,j)
      ENDDO
      xvenew(i) = xvenew(i)+parv(i)
    ENDIF
  ENDDO


! DEBUG: exit
! -----------
  IF (srfDebug) CALL cpu_exit(srName)

  RETURN

END SUBROUTINE itrf2etrf

END MODULE s_ITRF2ETRF
