! $Id: gpsleofix.f90 2022-10-28 09:51:28Z ychen $

!****s* fsi/gpsleofix.f90 *
!
! NAME
!   gpsleofix  -  satellites propagation to circular orbits
!
! SYNOPSIS
!    gpsleofix(nmax, r_gns, r_leo, r_coc, phase_sm, imon, r10, r20, teta, full_phase, iflag)
!
! DESCRIPTION
!
!   This subroutine calculates Doppler frequency shift (additional
!   to that in vacuum) associated with the increment of the zenith
!   angle (additional to that in vacuum) of the ray at one of the
!   satellites. Also some geometrical variables are computed.
!
!   Fix GPS and LEO radii to constants (at the point of occultation)
!   and correct the central angle and the phase between GPS and LEO
!   by continuation of rays, using externally calculated zenith angles.
!
! INPUTS
!    INTEGER,                    :: nmax       number of data points considered
!    REAL(wp), DIMENSION(:,:),   :: r_gns      GNSS coordinates (m)
!    REAL(wp), DIMENSION(:,:),   :: r_leo      LEO coordinates (m)
!    REAL(wp), DIMENSION(:),     :: r_coc      Centre curvature (m)
!    REAL(wp), DIMENSION(:,:),   :: phase_sm   excess phase after smooth (m)[ch, t]
!    INTEGER,                    :: imon       month (1-12)
!
! OUTPUT
!    REAL(wp),                   :: r10        absolute radius for GNSS (m)
!    REAL(wp),                   :: r20        absolute radius for LEO (m)
!    REAL(wp), DIMENSION(:),     :: teta       central angle between 1st and 2nd radius-vectors
!    REAL(wp), DIMENSION(:,:),   :: full_phase full phase (m) after propagating to circular orbit[ch, t]
!    INTEGER,                    :: iflag      0~teta increasing with time; 1~teta decreasing
!
! Reference: S.V.Sokolovskiy, Tracking tropospheric radio occulation signals
!            from low Earth orbit, Radio Science, Volume 36, Number 3, Pages
!            483-498, May/June 2001
!
! AUTHOR
!   S.V.Sokolovskiy, UCAR
!   update     2011/06/11  Feiqin Xie @ JPL (add parameter notation)
!   update     2014/03/10  Loknath Adhikari (update for airborne RO)
!   update     2019/02/28  Loknath Adhikari (modify for LEO)
!   Modified:  2022/10/28  Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE gpsleofix(nmax, r_gns, r_leo, r_coc, phase_sm, imon, r10, r20, teta, full_phase, iflag)

      USE typesizes, ONLY: wp => EightByteReal

      IMPLICIT NONE

      INTEGER,                  INTENT(in)    :: nmax      ! number of data points considered
      REAL(wp), DIMENSION(:,:), INTENT(in)    :: r_gns     ! GNSS coordinates (m)
      REAL(wp), DIMENSION(:,:), INTENT(in)    :: r_leo     ! LEO coordinates (m)
      REAL(wp), DIMENSION(:),   INTENT(in)    :: r_coc     ! Centre curvature (m)
      REAL(wp), DIMENSION(:,:), INTENT(in)    :: phase_sm  ! excess phase after smooth (m)[ch, t]
      INTEGER,                  INTENT(in)    :: imon      ! month
      REAL(wp),                 INTENT(out)   :: r10       ! absolute radius for GNSS (m)
      REAL(wp),                 INTENT(out)   :: r20       ! absolute radius for LEO (m)
      REAL(wp), DIMENSION(:),   INTENT(out)   :: teta      ! central angle between 1st and 2nd radius-vectors
      REAL(wp), DIMENSION(:,:), INTENT(out)   :: full_phase! full phase (m) after propagating to circular orbit[ch, t]
      INTEGER,                  INTENT(out)   :: iflag     ! 0~teta increasing with time; 1~teta decreasing


      INTEGER :: i, npoints, nxyz
      REAL(wp), dimension(size(r_leo,2))      :: r_gns_new ! Vect to GNSS from CoC
      REAL(wp), dimension(size(r_leo,2))      :: r_leo_new ! Vect to LEO from CoC
      REAL(wp) :: x1, y1, z1, x2, y2, z2, fi1, fi2, dfi2, glat, ras
      INTEGER  :: ierr
      REAL(wp) :: r10s, r20s, r1s, r2s, r1, r2, sp, ctet, dist
      REAL(wp) :: phi1, phi2, chi1, chi2, cf1, cf2, dph1, dph2, dtet1, dtet2

      npoints = SIZE(r_leo, 1)
      nxyz    = SIZE(r_leo, 2)

      iflag = 0
      r_gns_new = r_gns(nmax,:) - r_coc(:)
      r_leo_new = r_leo(nmax,:) - r_coc(:)

      r10s = SUM(r_gns_new**2)
      r20s = SUM(r_leo_new**2)
      r10  = SQRT(r10s)
      r20  = SQRT(r20s)

      DO i = 1, nmax
        x1 = r_gns(i,1) - r_coc(1)
        y1 = r_gns(i,2) - r_coc(2)
        z1 = r_gns(i,3) - r_coc(3)
        x2 = r_leo(i,1) - r_coc(1)
        y2 = r_leo(i,2) - r_coc(2)
        z2 = r_leo(i,3) - r_coc(3)

        ! compute azimuth angle
        CALL ray(imon, x1, y1, z1, x2, y2, z2, fi1, fi2, dfi2, glat, ras, ierr)

        r1s = x1**2 + y1**2 + z1**2
        r2s = x2**2 + y2**2 + z2**2
        r1 = SQRT(r1s)
        r2 = SQRT(r2s)
        sp = x1*x2 + y1*y2 + z1*z2
        ctet = sp / r1 / r2
        dist = SQRT(r1s + r2s - 2.0_wp * r1 * r2 * ctet)
        ! ---------------------------------------------------------------------
        ! For calculating phi1 and phi2 only
        ! ---------------------------------------------------------------------
        phi1 = ACOS( (r1s + dist**2 - r2s) / (2 * r1 * dist) )
        phi2 = ACOS( (r2s + dist**2 - r1s) / (2 * r2 * dist) )
        ! ---------------------------------------------------------------------
        chi1 = (r10s - r1s) / r1s
        chi2 = (r20s - r2s) / r2s
        cf1 = COS(fi1)
        cf2 = COS(fi2)
        dph1 = r1 * (chi1/2.0_wp/cf1 - chi1**2/8.0_wp/cf1**3 + chi1**3/16.0_wp/cf1**5)
        dph2 = r2 * (chi2/2.0_wp/cf2 - chi2**2/8.0_wp/cf2**3 + chi2**3/16.0_wp/cf2**5)
        dtet1 = ASIN(dph1 * SIN(fi1) / r10)
        dtet2 = ASIN(dph2 * SIN(fi2) / r20)
        teta(i) = ACOS(ctet) + dtet1 + dtet2
        full_phase(:, i) = dist + dph1 + dph2 + phase_sm(:, i)

        ! ---------------------------------------------------------------------
        !dnew = sqrt (r10s + r20s - 2 * r10 * r20 * dcos (teta (i) ) )
        !exphs = dist + dph1 + dph2 - dnew

      ENDDO

      DO i = 2, nmax
        IF (teta(i) <= teta(i - 1) ) iflag = 1
      ENDDO

END SUBROUTINE gpsleofix