      SUBROUTINE ELC(LAT0,LON0,LAT1,LON1,ROLL,PITCH,YAW,OFFSET,
     *                       LATC,LONC)
C
C___ parameter definition ___
C     IN
      REAL*8 LAT0(30),LON0(30),LAT1(30),LON1(30)
      REAL*8 ROLL(30),PITCH(30),YAW(30)
      INTEGER OFFSET
C     OUT
      REAL*8 LATC(30),LONC(30)
C     INTERNAL
      REAL*8 PI,RAD
      REAL*8 RE,F,EE,RN
      REAL*8 TI, TS, NP, SR, S0
      REAL*8 X1,Y1,Z1,UX1,UY1,UZ1,X2,Y2,Z2,UX2,UY2,UZ2
      REAL*8 X3,Y3,Z3,UX3,UY3,UZ3
      REAL*8 X4,Y4,Z4,UX4,UY4,UZ4
      REAL*8 UX5,UY5,UZ5,X6,Y6,Z6,UX6,UY6,UZ6,RATIO
      REAL*8 X7,Y7,Z7,UX7,UY7,UZ7,X8,Y8,Z8,UX8,UY8,UZ8
      REAL*8 X9,Y9,Z9,UX9,UY9,UZ9
      REAL*8 UX10,UY10,UZ10,X11,Y11,Z11,UX11,UY11,UZ11
      REAL*8 X12,Y12,Z12,UX12,UY12,UZ12
      REAL*8 UXN,UYN,UZN
      REAL*8 ANGLE12,ANGLE67,ANGLE510
      REAL*8 UXYAW,UYYAW,UZYAW
      REAL*8 LATP,LONP,XP,YP,ZP,UXP,UYP,UZP
      REAL*8 XROLL,YROLL,ZROLL,UXROLL,UYROLL,UZROLL
      REAL*8 UXPITCH,UYPITCH,UZPITCH
      REAL*8 ANGLENP,XO,YO,ZO,altitude,SCAN
      REAL*8 CM(3,3)
      REAL*8 DR,DP,DY
      REAL*8 UXCPP,UYCPP,UZCPP,UXCP,UYCP,UZCP
      REAL*8 KA,KB,KC,K1,K2,K,XPC,YPC,ZPC
      INTEGER P
C    PI
      PI = 3.14159265359D0
      RAD = PI/180.D0
C    Earth shape parameters
      RE = 6378.135D0
      F = 1.D0 / 298.25D0
      EE = 2.D0 * F - F * F
C    Instrument parameters
C    DuRATIOn between adjacent IFOVs (sec)
      TI = 0.2025D0
C    Time interval between adjacent SCANs (sec)
      TS = 8D0
C    SCAN position of nadir
      PC = 15.5D0
C    SCAN angle interval between adjacent positions (degree)
      SR = (3.D0 + 1.D0 / 3.D0)
C    SCAN angle of position #0 (degree)
      S0 = SR * PC
C
C___ correct earth location data ___
C
C    ___ scan n ___
C    ___ vector 1 ___
      RN = RE / DSQRT(1.D0-EE*DSIN(LAT0(15)*RAD)*DSIN(LAT0(15)*RAD))
      X1 = RN * DCOS(LAT0(15)*RAD) * DCOS((360.D0+LON0(15))*RAD)
      Y1 = RN * DCOS(LAT0(15)*RAD) * DSIN((360.D0+LON0(15))*RAD)
      Z1 = RN * (1.D0-EE) * DSIN(LAT0(15)*RAD)
      UX1 = X1 / DSQRT(X1 * X1 + Y1 * Y1 + Z1 * Z1)
      UY1 = Y1 / DSQRT(X1 * X1 + Y1 * Y1 + Z1 * Z1)
      UZ1 = Z1 / DSQRT(X1 * X1 + Y1 * Y1 + Z1 * Z1)
C    ___ vector 2 ___
      RN = RE / DSQRT(1.D0-EE*DSIN(LAT0(16)*RAD)*DSIN(LAT0(16)*RAD))
      X2 = RN * DCOS(LAT0(16)*RAD) * DCOS(LON0(16)*RAD)
      Y2 = RN * DCOS(LAT0(16)*RAD) * DSIN(LON0(16)*RAD)
      Z2 = RN * (1.D0-EE) * DSIN(LAT0(16)*RAD)
      UX2 = X2 / DSQRT(X2 * X2 + Y2 * Y2 + Z2 * Z2)
      UY2 = Y2 / DSQRT(X2 * X2 + Y2 * Y2 + Z2 * Z2)
      UZ2 = Z2 / DSQRT(X2 * X2 + Y2 * Y2 + Z2 * Z2)
C    ___ vector 3 ___
      X3 = UY1 * UZ2 - UZ1 * UY2
      Y3 = UZ1 * UX2 - UX1 * UZ2
      Z3 = UX1 * UY2 - UY1 * UX2
      UX3 = X3 / DSQRT(X3 * X3 + Y3 * Y3 + Z3 * Z3)
      UY3 = Y3 / DSQRT(X3 * X3 + Y3 * Y3 + Z3 * Z3)
      UZ3 = Z3 / DSQRT(X3 * X3 + Y3 * Y3 + Z3 * Z3)
C    ___ vector 4 ___
      X4 = UY3 * UZ1 - UZ3 * UY1
      Y4 = UZ3 * UX1 - UX3 * UZ1
      Z4 = UX3 * UY1 - UY3 * UX1
      UX4 = X4 / DSQRT(X4 * X4 + Y4 * Y4 + Z4 * Z4)
      UY4 = Y4 / DSQRT(X4 * X4 + Y4 * Y4 + Z4 * Z4)
      UZ4 = Z4 / DSQRT(X4 * X4 + Y4 * Y4 + Z4 * Z4)
C    ___ angle 1 to 2 ___
      ANGLE12 = DACOS(UX1 * UX2 + UY1 * UY2 + UZ1 * UZ2)
C    ___ vector 5 : nadir ___
      UX5 = UX1 * DCOS(ANGLE12 * 0.5D0) + UX4 * DSIN(ANGLE12 * 0.5D0)
      UY5 = UY1 * DCOS(ANGLE12 * 0.5D0) + UY4 * DSIN(ANGLE12 * 0.5D0)
      UZ5 = UZ1 * DCOS(ANGLE12 * 0.5D0) + UZ4 * DSIN(ANGLE12 * 0.5D0)
C
C    ___ scan n+1 ___
C    ___ vector 6 ___
      RN = RE / DSQRT(1.D0-EE*DSIN(LAT1(15)*RAD)*DSIN(LAT1(15)*RAD))
      X6 = RN * DCOS(LAT1(15)*RAD) * DCOS(LON1(15)*RAD)
      Y6 = RN * DCOS(LAT1(15)*RAD) * DSIN(LON1(15)*RAD)
      Z6 = RN * (1.D0-EE) * DSIN(LAT1(15)*RAD)
      UX6 = X6 / DSQRT(X6 * X6 + Y6 * Y6 + Z6 * Z6)
      UY6 = Y6 / DSQRT(X6 * X6 + Y6 * Y6 + Z6 * Z6)
      UZ6 = Z6 / DSQRT(X6 * X6 + Y6 * Y6 + Z6 * Z6)
C    ___ vector 7 ___
      RN = RE / DSQRT(1-EE*DSIN(LAT1(16)*RAD)*DSIN(LAT1(16)*RAD))
      X7 = RN * DCOS(LAT1(16)*RAD) * DCOS(LON1(16)*RAD)
      Y7 = RN * DCOS(LAT1(16)*RAD) * DSIN(LON1(16)*RAD)
      Z7 = RN * (1.D0-EE) * DSIN(LAT1(16)*RAD)
      UX7 = X7 / DSQRT(X7 * X7 + Y7 * Y7 + Z7 * Z7)
      UY7 = Y7 / DSQRT(X7 * X7 + Y7 * Y7 + Z7 * Z7)
      UZ7 = Z7 / DSQRT(X7 * X7 + Y7 * Y7 + Z7 * Z7)
C    ___ vector 8 ___
      X8 = UY6 * UZ7 - UZ6 * UY7
      Y8 = UZ6 * UX7 - UX6 * UZ7
      Z8 = UX6 * UY7 - UY6 * UX7
      UX8 = X8 / DSQRT(X8 * X8 + Y8 * Y8 + Z8 * Z8)
      UY8 = Y8 / DSQRT(X8 * X8 + Y8 * Y8 + Z8 * Z8)
      UZ8 = Z8 / DSQRT(X8 * X8 + Y8 * Y8 + Z8 * Z8)
C    ___ vector 9 ___
      X9 = UY8 * UZ6 - UZ8 * UY6
      Y9 = UZ8 * UX6 - UX8 * UZ6
      Z9 = UX8 * UY6 - UY8 * UX6
      UX9 = X9 / DSQRT(X9 * X9 + Y9 * Y9 + Z9 * Z9)
      UY9 = Y9 / DSQRT(X9 * X9 + Y9 * Y9 + Z9 * Z9)
      UZ9 = Z9 / DSQRT(X9 * X9 + Y9 * Y9 + Z9 * Z9)
C    ___ angle 6 to 7 ___
      ANGLE67 = DACOS(UX6 * UX7 + UY6 * UY7 + UZ6 * UZ7)
C    ___ vector 10 : nadir ___
      UX10 = UX6 * DCOS(ANGLE67 * 0.5D0) + UX9 * DSIN(ANGLE67 * 0.5D0)
      UY10 = UY6 * DCOS(ANGLE67 * 0.5D0) + UY9 * DSIN(ANGLE67 * 0.5D0)
      UZ10 = UZ6 * DCOS(ANGLE67 * 0.5D0) + UZ9 * DSIN(ANGLE67 * 0.5D0)
C
C    ___ vector 11 ___
      X11 = UY5 * UZ10 - UZ5 * UY10
      Y11 = UZ5 * UX10 - UX5 * UZ10
      Z11 = UX5 * UY10 - UY5 * UX10
      UX11 = X11 / DSQRT(X11 * X11 + Y11 * Y11 + Z11 * Z11)
      UY11 = Y11 / DSQRT(X11 * X11 + Y11 * Y11 + Z11 * Z11)
      UZ11 = Z11 / DSQRT(X11 * X11 + Y11 * Y11 + Z11 * Z11)
C    ___ vector 12 ___
      X12 = UY11 * UZ5 - UZ11 * UY5
      Y12 = UZ11 * UX5 - UX11 * UZ5
      Z12 = UX11 * UY5 - UY11 * UX5
      UX12 = X12 / DSQRT(X12 * X12 + Y12 * Y12 + Z12 * Z12)
      UY12 = Y12 / DSQRT(X12 * X12 + Y12 * Y12 + Z12 * Z12)
      UZ12 = Z12 / DSQRT(X12 * X12 + Y12 * Y12 + Z12 * Z12)
C    ___ angle 5 to 10 ___
      ANGLE510 = DACOS(UX5 * UX10 + UY5 * UY10 + UZ5 * UZ10)
C    ___ correction for each position ___
      DO 1000 P = 1, 30
C        ___ nadir ___
        RATIO = (DBLE(P) - PC) * TI / TS + DBLE(OFFSET)
        UXN = UX5*DCOS(ANGLE510 * RATIO) + UX12*DSIN(ANGLE510 * RATIO)
        UYN = UY5*DCOS(ANGLE510 * RATIO) + UY12*DSIN(ANGLE510 * RATIO)
        UZN = UZ5*DCOS(ANGLE510 * RATIO) + UZ12*DSIN(ANGLE510 * RATIO)
C        ___ yaw (+XP) vector ___
        UXYAW = -UXN
        UYYAW = -UYN
        UZYAW = -UZN
C        ___ scene vector ___
        SCAN = (-SR * DBLE(P) + S0) * RAD
       IF(OFFSET.EQ.0) THEN
          RN = RE / DSQRT(1.D0-EE*DSIN(LAT0(P)*RAD)*DSIN(LAT0(P)*RAD))
          XP = RN * DCOS(LAT0(P)*RAD) * DCOS(LON0(P)*RAD)
          YP = RN * DCOS(LAT0(P)*RAD) * DSIN(LON0(P)*RAD)
          ZP = RN * (1.D0-EE) * DSIN(LAT0(P)*RAD)
       ELSEIF(OFFSET.EQ.1) THEN
          RN = RE / DSQRT(1.D0-EE*DSIN(LAT1(P)*RAD)*DSIN(LAT1(P)*RAD))
          XP = RN * DCOS(LAT1(P)*RAD) * DCOS(LON1(P)*RAD)
          YP = RN * DCOS(LAT1(P)*RAD) * DSIN(LON1(P)*RAD)
          ZP = RN * (1.D0-EE) * DSIN(LAT1(P)*RAD)
       ENDIF
        UXP = XP / DSQRT(XP * XP + YP * YP + ZP * ZP)
        UYP = YP / DSQRT(XP * XP + YP * YP + ZP * ZP)
        UZP = ZP / DSQRT(XP * XP + YP * YP + ZP * ZP)
C        ___ roll (+YP) vector ___
        XROLL = (UYYAW * UZP - UZYAW * UYP) *(-DSIGN(1.D0,SCAN))
        YROLL = (UZYAW * UXP - UXYAW * UZP) *(-DSIGN(1.D0,SCAN))
        ZROLL = (UXYAW * UYP - UYYAW * UXP) *(-DSIGN(1.D0,SCAN))
        UXROLL = XROLL/DSQRT(XROLL*XROLL+YROLL*YROLL+ZROLL*ZROLL)
        UYROLL = YROLL/DSQRT(XROLL*XROLL+YROLL*YROLL+ZROLL*ZROLL)
        UZROLL = ZROLL/DSQRT(XROLL*XROLL+YROLL*YROLL+ZROLL*ZROLL)
C        ___ pitch (+ZP) vector ___
        UXPITCH = UYYAW * UZROLL - UZYAW * UYROLL
        UYPITCH = UZYAW * UXROLL - UXYAW * UZROLL
        UZPITCH = UXYAW * UYROLL - UYYAW * UXROLL
C        ___ satellite point ___
        ANGLENP = DACOS(UXN * UXP + UYN * UYP + UZN * UZP)
        XO =(DSIN(ANGLENP)/DTAN(DABS(SCAN))+DCOS(ANGLENP))*
     .       DSQRT(XP*XP+YP*YP+ZP*ZP)*UXN
        YO =(DSIN(ANGLENP)/DTAN(DABS(SCAN))+DCOS(ANGLENP))*
     .       DSQRT(XP*XP+YP*YP+ZP*ZP)*UYN
        ZO =(DSIN(ANGLENP)/DTAN(DABS(SCAN))+DCOS(ANGLENP))*
     .       DSQRT(XP*XP+YP*YP+ZP*ZP)*UZN
        ALTITUDE = DSQRT(XO * XO + YO * YO + ZO * ZO) - RE
C        ___ corrected pointing vector ___
        DR = DBLE(ROLL(P))
        DP = DBLE(PITCH(P))
        DY = DBLE(YAW(P))

        CM(1, 1) = DCOS(DP)*DCOS(DR)
        CM(1, 2) = DCOS(DP)*DSIN(DR)*DSIN(DY)-DSIN(DP)*DCOS(DY)
        CM(1, 3) = DCOS(DP)*DSIN(DR)*DCOS(DY)+DSIN(DP)*DSIN(DY)
        CM(2, 1) = DSIN(DP)*DCOS(DR)
        CM(2, 2) = DSIN(DP)*DSIN(DR)*DSIN(DY)+DCOS(DP)*DCOS(DY)
        CM(2, 3) = DSIN(DP)*DSIN(DR)*DCOS(DY)-DCOS(DP)*DSIN(DY)
        CM(3, 1) = -DSIN(DR)
        CM(3, 2) = DCOS(DR) * DSIN(DY)
        CM(3, 3) = DCOS(DR) * DCOS(DY)
       UXCPP = CM(1, 1) * DCOS(SCAN) + CM(1, 3) * DSIN(SCAN)
        UYCPP = CM(2, 1) * DCOS(SCAN) + CM(2, 3) * DSIN(SCAN)
        UZCPP = CM(3, 1) * DCOS(SCAN) + CM(3, 3) * DSIN(SCAN)
        UXCP = UXYAW * UXCPP + UXROLL * UYCPP + UXPITCH * UZCPP
        UYCP = UYYAW * UXCPP + UYROLL * UYCPP + UYPITCH * UZCPP
        UZCP = UZYAW * UXCPP + UZROLL * UYCPP + UZPITCH * UZCPP
C        ___ corrected view vector ___
        KA = (1-F)*(1-F)*(UXCP*UXCP+UYCP*UYCP)+UZCP*UZCP
        KB = (1-F)*(1-F)*(XO*UXCP+YO*UYCP)+ZO*UZCP
        KC = (1-F)*(1-F)*(XO*XO+YO*YO-RE*RE)+ZO*ZO
        K1 = (-KB+DSQRT(KB*KB-KA*KC))/KA
        K2 = (-KB-DSQRT(KB*KB-KA*KC))/KA
        If(DABS(K1).GT.DABS(K2)) Then
            K = K2
        Else
            K = K1
        End If
C        ___ corrected scene point ___
        XPC = XO + K * UXCP
        YPC = YO + K * UYCP
        ZPC = ZO + K * UZCP
        LATC(P) = DATAN(ZPC/((1-F)*(1-F)*DSQRT(XPC*XPC+YPC*YPC)))
        If(XPC.NE.0) Then
            LONC(P) = DATAN(YPC/XPC)
            If(XPC.LT.0.AND.YPC.GE.0) LONC(P) = LONC(P) + PI
            If(XPC.LT.0.AND.YPC.LT.0) LONC(P) = LONC(P) - PI
        Else
            If(YPC.GT.0) Then
                LONC(P) = PI / 2
            Else
                LONC(P) = -PI / 2
            End If
        End If
 1000  CONTINUE
C
      RETURN
      END
