      SUBROUTINE ORBINT(IFLAG,T,DLAT,DLONG,HT)
C     SUBROUTINE ORBINT: INPUT 4 WORLD MAP POINTS (GMT,LAT,LONG,HT)
C     NEAREST THE IONOGRAM START TIME. INTERPOLATE TO FIND SATELLITE
C     POSITION AT TIME T, WHICH MAY BE THE START TIME OR THE TIME OF A
C     SCALED FREQUENCY (WHEN SATELLITE POSITION IS ALLOWED TO VARY OVER
C     THE IONOGRAM). INTERPOLATION USES AITKEN S THIRD-ORDER REPEATED
C     LINEAR INTERPOLATION METHOD (CONTAINED IN THE FUNCTION TAB).
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*4 DLAT,DLONG,HT
      DIMENSION GMT(4),DLATL(4),DLONGL(4),HGHT(4),XP(4),YP(4)
      SAVE
      DATA RAD/57.29577951D0/
      IF (IFLAG .NE. 0) GO TO 100
      READ (3,1) (GMT(K),DLATL(K),DLONGL(K),HGHT(K),K=1,4)
C     IFLAG=0 FOR INITIAL CALL PER IONOGRAM TO ORBINT. CONVERT ALL TIMES
C     TO DECIMAL.  IF NEEDED, CONVERT LONGITUDES TO THE RANGE 0-360, AND
C     FORM ARRAYS XP AND YP (POLAR PROJECTIONS OF DLATL,DLONGL) FOR
C     LATITUDE AND LONGITUDE INTERPOLATION WHEN LATITUDE DOES NOT CHANGE
C     SIGN IN THE TIME INTERVAL GMT(1) TO GMT(4).
      ILONG=0
      IF(DABS(DLONGL(1)).LT.90.D0) ILONG=1
      IDIF=0
      DO 400 J=2,4
      IF (GMT(J) .LT. GMT(J-1)) GMT(J)=GMT(J)+240000.D0
      IF (DLATL(J)*DLATL(J-1) .LT. 0.0D0) IDIF=1
  400 CONTINUE
      DO 200 K=1,4
      XX=GMT(K)
      GMT(K)=0.0D0
      DO 300 J=1,3
      XX=XX/100.0D0
      JJ=XX
      GMT(K)=GMT(K)+(XX-JJ)*100.0D0/(60.0D0**(3-J))
  300 XX=JJ
      IF((IDIF*ILONG).EQ.1) GO TO 200
      IF (DLONGL(K) .LT. 0.0D0) DLONGL(K)=360.0D0+DLONGL(K)
      XX=90.0D0-DABS(DLATL(K))
      YY=DLONGL(K)/RAD
      XP(K)=XX*DCOS(YY)
      YP(K)=XX*DSIN(YY)
  200 CONTINUE
      SIGN=1.0D0
      IF (DLATL(1) .LT. 0.0D0) SIGN=-1.0D0
  100 IF (T .LT. GMT(1)) T=T+24.D0
      IF (IDIF .EQ. 1) GO TO 500
C     LATITUDE DOES NOT CHANGE SIGN IN THIS TIME INTERVAL. INTERPOLATE
C     ON XP AND YP TO GET LATITUDE AND LONGITUDE AT TIME T.
      XX=TAB(GMT,XP,4,3,T)
      YY=TAB(GMT,YP,4,3,T)
      DLAT=(90.0D0-DSQRT(XX*XX+YY*YY))*SIGN
      DLONG=DATAN2(YY,XX)*RAD
      GO TO 600
C     LATITUDE CHANGES SIGN IN THIS TIME INTERVAL.
  500 DLAT=TAB(GMT,DLATL,4,3,T)
      DLONG=TAB(GMT,DLONGL,4,3,T)
C     CONVERT LONGITUDE TO RANGE -180 TO +180.   COMPUTE SATELLITE HT.
  600 IF(DLONG.GT.180.) DLONG=DLONG-360.
      HT=TAB(GMT,HGHT,4,3,T)
      RETURN
    1 FORMAT (8F10.4)
      END
