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

MODULE s_ITRF2ITRF

  IMPLICIT NONE

#include "COM_CPU.inc"

CONTAINS

SUBROUTINE itrf2itrf(timcrd,datin,datout,xstat,xstnew,xvel,xvenew)

! -------------------------------------------------------------------------
! Purpose:    Transformation from ITRFyy 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_
!             22-Dec-2020 RD: Transform D_CONST to a protected structure CONST
!             28-Nov-2023 SL: ITRF88 and ITRF2014 w.r.t. ITRF2020 added
!             28-Nov-2023 SL: Print transformation parameters
!             28-Nov-2023 SL: Support of velocity transformation
!             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=*)                   :: datin      ! Source datum string
  CHARACTER(LEN=*)                   :: datout     ! Target datum string
  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 = 'itrf2itrf'
  INTEGER(i4b),     PARAMETER        :: ntrf = 14  ! Nbr of supported ITRFs

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

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

  INTEGER(i4b)                       :: i, j

  LOGICAL,                      SAVE :: first = .TRUE.


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


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

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


! Transformation parameters from ITRF2020 to past ITRFs
! -----------------------------------------------------
!                  Tx[mm]   Ty[mm]   Tz[mm]  D[ppb]   Rx[mas]   Ry[mas]   Rz[mas]   t0[y]
! ITRF2020
  param(14,:) = (/  0.0d0,   0.0d0,   0.0d0,  0.00d0,  0.000d0,  0.000d0,  0.000d0, 2015.0d0/)
  rates(14,:) = (/  0.0d0,   0.0d0,   0.0d0,  0.00d0,  0.000d0,  0.000d0,  0.000d0/)
! ITRF2014
  param(13,:) = (/ -1.4d0,  -0.9d0,   1.4d0, -0.42d0,  0.000d0,  0.000d0,  0.000d0, 2015.0d0/)
  rates(13,:) = (/  0.0d0,  -0.1d0,   0.2d0,  0.00d0,  0.000d0,  0.000d0,  0.000d0/)
! ITRF2008
  param(12,:) = (/  0.2d0,   1.0d0,   3.3d0, -0.29d0,  0.000d0,  0.000d0,  0.000d0, 2015.0d0/)
  rates(12,:) = (/  0.0d0,  -0.1d0,   0.1d0,  0.03d0,  0.000d0,  0.000d0,  0.000d0/)
! ITRF2005
  param(11,:) = (/  2.7d0,   0.1d0,  -1.4d0,  0.65d0,  0.000d0,  0.000d0,  0.000d0, 2015.0d0/)
  rates(11,:) = (/  0.3d0,  -0.1d0,   0.1d0,  0.03d0,  0.000d0,  0.000d0,  0.000d0/)
! ITRF2000
  param(10,:) = (/ -0.2d0,   0.8d0, -34.2d0,  2.25d0,  0.000d0,  0.000d0,  0.000d0, 2015.0d0/)
  rates(10,:) = (/  0.1d0,   0.0d0,  -1.7d0,  0.11d0,  0.000d0,  0.000d0,  0.000d0/)
! ITRF97
  param( 9,:) = (/  6.5d0,  -3.9d0, -77.9d0,  3.98d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 9,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF96
  param( 8,:) = (/  6.5d0,  -3.9d0, -77.9d0,  3.98d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 8,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF94
  param( 7,:) = (/  6.5d0,  -3.9d0, -77.9d0,  3.98d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 7,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF93
  param( 6,:) = (/-65.8d0,   1.9d0, -71.3d0,  4.47d0, -3.360d0, -4.330d0,  0.750d0, 2015.0d0/)
  rates( 6,:) = (/ -2.8d0,  -0.2d0,  -2.3d0,  0.12d0, -0.110d0, -0.190d0,  0.070d0/)
! ITRF92
  param( 5,:) = (/ 14.5d0,  -1.9d0, -85.9d0,  3.27d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 5,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF91
  param( 4,:) = (/ 26.5d0,  12.1d0, -91.9d0,  4.67d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 4,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF90
  param( 3,:) = (/ 24.5d0,   8.1d0,-107.9d0,  4.97d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 3,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF89
  param( 2,:) = (/ 29.5d0,  32.1d0,-145.9d0,  8.37d0,  0.000d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 2,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)
! ITRF88
  param( 1,:) = (/ 24.5d0,  -3.9d0,-169.9d0, 11.47d0,  0.100d0,  0.000d0,  0.360d0, 2015.0d0/)
  rates( 1,:) = (/  0.1d0,  -0.6d0,  -3.1d0,  0.12d0,  0.000d0,  0.000d0,  0.020d0/)


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

  DO j=1,ntrf
    IF(datin(len_trim(datin)-3:len_trim(datin)) == 'R3_0') THEN
      IF(INDEX(datin(len_trim(datin)-5:len_trim(datin)-4),trflst(j))==0) CYCLE
    ELSEIF(datin(len_trim(datin)-1:len_trim(datin)) == 'R3') THEN
      IF(INDEX(datin(len_trim(datin)-3:len_trim(datin)-2),trflst(j))==0) CYCLE
    ELSEIF(datin(len_trim(datin):len_trim(datin)) == 'P') THEN
      IF(INDEX(datin(len_trim(datin)-2:len_trim(datin)-1),trflst(j))==0) CYCLE
    ELSEIF(datin(len_trim(datin)-1:len_trim(datin)) == '_0') THEN
      IF(INDEX(datin(len_trim(datin)-3:len_trim(datin)-2),trflst(j))==0) CYCLE
    ELSE
      IF(INDEX(datin(len_trim(datin)-1:len_trim(datin)),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)') 'Input datum: ',trim(datin)
    CALL exitRc(2)
  ENDIF

  i = 0

  DO j=1,ntrf
    IF(INDEX(datout(len_trim(datout)-1:len_trim(datout)),trflst(j))==0) CYCLE
    xx   = trflst(j)
    paro = param(j,:)
    rato = rates(j,:)
    doy  = (param(j,8)-REAL(INT(param(j,8))))*365.25d0+1d0
    dto  = (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)') 'Output datum: ',trim(datout)
    CALL exitRc(2)
  ENDIF

  IF(ABS(pari(8)-paro(8))>0d0) THEN
    WRITE(lfnErr,"(/,' *** SR ',A,':',A)") srname, &
         ' Difference in the reference epoch found'
    WRITE(lfnErr,'(8X,A,F7.2)') 'Reference epoch input frame : ',pari(8)
    WRITE(lfnErr,'(8X,A,F7.2)') 'Reference epoch output frame: ',paro(8)
    CALL exitRc(2)
  ENDIF

  IF(first) THEN
    WRITE(lfnPrt,*)
    WRITE(lfnPrt,'(1X,A,A,A,A,A)') &
         'ITRF-',yy,' to ITRF-',xx,':'
    WRITE(lfnPrt,'(1X,18("-"))')
    WRITE(lfnPrt,*)

    WRITE(lfnPrt,'(1X,73("-"))')
    WRITE(lfnPrt,'(1X,3A8,A9,3A10,2X,A8)') &
         'dTx','dTy','dTz', &
         'dD', &
         'dRx','dRy','dRz', &
         '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)') &
         paro(1)-pari(1),paro(2)-pari(2),paro(3)-pari(3), &
         paro(4)-pari(4), &
         paro(5)-pari(5),paro(6)-pari(6),paro(7)-pari(7), &
         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)') &
         rato(1)-rati(1),rato(2)-rati(2),rato(3)-rati(3), &
         rato(4)-rati(4), &
         rato(5)-rati(5),rato(6)-rati(6),rato(7)-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) = ((paro(i)+rato(i)*dto)-(pari(i)+rati(i)*dti))*1d-3
      parv(i) = (rato(i)-rati(i))*1d-3
    ELSEIF(i==4) THEN
      parc(i) = ((paro(i)+rato(i)*dto)-(pari(i)+rati(i)*dti))*1d-9
      parv(i) = (rato(i)-rati(i))*1d-9
    ELSEIF(i==5.OR.i==6.OR.i==7) THEN
      parc(i) = ((paro(i)+rato(i)*dto)-(pari(i)+rati(i)*dti))*1d-3/CONST%ars
      parv(i) = (rato(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 itrf2itrf

END MODULE s_ITRF2ITRF
