OceanToRayX Function

FUNCTION OceanToRayX( x, xs, tradial, t, snapDim )
!! Transform ocean coordinates to ray coordinates

  ! LP: Going back and forth through the coordinate transform won't
  ! always keep the precise value, so we may have to finesse the floats.
  ! Valid values of snapDim:
  ! -2: Snap to X or Y unspecified
  ! -1: No snap
  !  0: Snap to X
  !  1: Snap to Y
  !  2: Snap to Z

  REAL ( KIND=8 ) :: OceanToRayX( 2 )
  REAL ( KIND=8 ), INTENT( IN ) :: x( 3 ), xs( 3 ), tradial( 2 ), t( 2 )
  INTEGER, INTENT( IN ) :: snapDim
  REAL ( KIND=8 ) :: ret( 2 ), x_back( 3 ), wantdir( 2 ), errdir( 2 )
  LOGICAL :: correctdir( 2 )
  INTEGER :: i

  ! Depth always transfers perfectly--not changed.
  ret( 2 ) = x( 3 )
  ! For range, use larger dimension--this avoids divide-by-zero or divide
  ! by a small number causing accuracy problems.
  IF ( ABS( tradial( 1 ) ) >= ABS( tradial( 2 ) ) ) THEN
     ret( 1 ) = ( x( 1 ) - xs( 1 ) ) / tradial( 1 )
  ELSE
     ret( 1 ) = ( x( 2 ) - xs( 2 ) ) / tradial( 2 )
  END IF
  IF ( snapDim < -2 .OR. snapDim == -1 .OR. snapDim >= 2 ) THEN
     ! Either:
     ! snapDim out of range (won't happen, but want to help compiler)
     ! No snap selected--this is the best estimate
     ! Snap to Z--Z already perfect, this is the best estimate
     OceanToRayX = ret
     RETURN
  END IF
  ! Only do this iteration a few times, then give up.
  DO i = 1, 4
     ! Go back from 2D to 3D, compare to original x.
     x_back = RayToOceanX( ret, xs, tradial )
     ! If we can't be on the boundary, we want to be slightly forward of
     ! the boundary, measured in terms of the ray tangent range. This also
     ! encompasses cases where one component exactly matches (errdir.x_or_y
     ! == RL(0.0)). For both of these values, only the sign matters.
     wantdir = tradial * t( 1 )
     errdir = x_back( 1 : 2 ) - x( 1 : 2 )
     correctdir = ( wantdir * errdir ) >= 0.0D0
     IF ( ( snapDim == 0 .AND. correctdir( 1 ) ) .OR. &
          ( snapDim == 1 .AND. correctdir( 2 ) ) .OR. &
          ( correctdir( 1 ) .AND. correctdir( 2 ) ) ) THEN
        OceanToRayX = ret
        RETURN
     END IF
     ! Move to the next floating point value for ret, in the direction
     ! of the ray tangent. How do we know this will actually produce a
     ! new value for x_back? If the scale of the floating point steps
     ! around x_back is larger than that of ret (several values of ret
     ! map to the same value of x_back after the transform), there
     ! should be no problem in finding a value that maps exactly, so
     ! we would not get here.
     ret( 1 ) = NEAREST( ret( 1 ), t( 1 ) )
  END DO
  WRITE( PRTFile, * ) 'Warning in OceanToRayX: Failed to transform 3D -> 2D -> 3D consistently'
  OceanToRayX = ret
END FUNCTION OceanToRayX