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

!****s* fsi/ray.f90 *
!
! NAME
!   ray  -  This subroutine calculates parameters of ray
!
! SYNOPSIS
!    call ray(imonth, x1, y1, z1, x2, y2, z2, fi1, fi2, dfi2, glat, ras, ierr)
!
! DESCRIPTION
!    This subroutine calculates parameters of ray (zenith angles,
!    altitude of ray asympthote over the reference ellipsoid, and
!    latitude of the tangent point) for given positions of transmitter
!    and receiver and bending angle climatology.
!
! INPUTS
!    INTEGER,   :: imonth                    month (1-12)
!    REAL(wp),  :: x1, y1, z1, x2, y2, z2    x1,y1,z1 - GPS position (km)
!                                            x2,y2,z2 - LEO position (km)
!
! OUTPUT
!    REAL(wp),  :: fi1, fi2, dfi2, glat, ras
!                    fi1      - zenith angles of the ray at GPS
!                    fi2      - zenith angles of the ray at LEO
!                    dfi2     - increment of arrival angle of the ray at 2nd satellite
!                    glat     - latitude of the tangent point (deg)
!                    ras      - altitude of ray asympthote over the reference ellipsoid
!    INTEGER,   :: ierr   flag 0 = OK, 1 = error
!
! 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
!   Modified: 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 ray(imonth, x1, y1, z1, x2, y2, z2, fi1, fi2, dfi2, glat, ras, ierr)

      USE typesizes, ONLY: wp => EightByteReal
      USE ropp_utils, ONLY: vector_product, vector_angle, rotate
      USE ropp_pp_constants, ONLY: Pi

      IMPLICIT NONE
      INTEGER,  INTENT(in)    :: imonth
      REAL(wp), INTENT(in)    :: x1, y1, z1, x2, y2, z2
      REAL(wp), INTENT(out)   :: fi1, fi2, dfi2, glat, ras
      INTEGER,  INTENT(out)   :: ierr

      REAL(wp), DIMENSION(3)  :: v1(3), v2(3), vs(3), vt(3)
      REAL(wp) :: r1, r2, omega, talpha, alpha, a0, teta, tet, dtet, der, dt, glat1
      REAL(wp) :: a, b, b1, amin, bmax
      INTEGER  :: ic, iflag1, iflag2

      v1(1) = x1
      v1(2) = y1
      v1(3) = z1
      v2(1) = x2
      v2(2) = y2
      v2(3) = z2

      vs = vector_product(v2, v1)

      r1 = SQRT(SUM(v1(:)**2))
      r2 = SQRT(SUM(v2(:)**2))

      teta = vector_angle(v2, v1)
      omega = Pi - teta
      talpha = r1*Sin(omega) / (r2 + r1*Cos(omega))
      a0 = r2*talpha / Sqrt(1.0_wp + talpha**2)

      fi2 = dasin (a0 / r2)

      vt = rotate(v2, vs, pi / 2.0_wp - fi2)
      glat = 90.0_wp - 180.0_wp/Pi * acos(vt(3)/sqrt(sum(vt(:)**2)))

      iflag1 = 0
      iflag2 = 0
      ierr = 0

      a = a0

      ic = 0

      DO WHILE (ic < 1000 )

         CALL benmod(imonth, glat, a, b, b1, amin, bmax)
         dtet = pi - dasin(a/r1) - dasin(a/r2) + b - teta
         IF (dabs(dtet) > 1.0d-5) THEN
           der = - 1.0_wp / SQRT(r1**2 - a**2) - 1.0_wp/ SQRT(r2**2 - a**2) + b1
           a = a - dtet / der
         ELSE
           iflag1 = 1
         ENDIF

         IF (a < amin) THEN
           a   = amin
           b   = bmax
           tet = Pi-ASIN(a/r1)-ASIN(a/r2) + b
           dt  = teta-tet
         ELSE
           dt = 0.d0
         ENDIF

         fi1  = ASIN(a/r1)
         fi2  = ASIN(a/r2)

         dfi2 = fi2 - ASIN(a0/r2)

         alpha = Pi/2.0_wp-fi2 + b/2.0_wp + dt/2.0_wp

         vt = rotate(v2, vs, alpha)
         glat1 = 90.0_wp - 180.0_wp/Pi * acos(vt(3)/sqrt(sum(vt(:)**2)))
         IF (dabs(glat1-glat) > 1.0d-3) THEN
           glat=glat1
         ELSE
           iflag2 = 1
         ENDIF
         ic=ic+1

         ! Both conditions are met
         IF (iflag1 > 0 .AND. iflag2 > 0 ) THEN
           exit
         ENDIF

      ENDDO

      IF (ic > 1000) THEN
         ierr = 1
      ENDIF

      ras = a - amin

END SUBROUTINE ray