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))*DSIN(LAT1(16))) 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