|
|
NOAA GVI GUIDEAppendix L |
Introduction Page,
GVI Guide TOC, Acronyms
This appendix contains the FORTRAN source code for a subroutine (RELTAZ) which calculates the relative azimuth angle. Two parameters, TIME (time after satellite launch) and DT (equatorial crossing) must be calculated prior to calling subroutine RELTAZ.
To calculate TIME which is a function of a specific satellite, use the following equations:
where, YR is the four digit qualifier for the year and MON is the numerical representation for month (e.g., January would be represented by "1", February would be "2", etc.).
To calculate the equatorial crossing, DT, which is also a function of satellite and node, use the following polynomial:
Note that DT is in hours on a 24-hour clock. The equations given above are for the afternoon crossing. If a user wants to use PM hours on a 12-hour clock, then 12 should be subtracted from the first term of each equation.
Once the other input variables are given (pixel latitude, pixel longitude, the negative of the scan angle, the Julian day, and the actual solar zenith angle), then the function RELTAZ can be called. See the actual source code for details about allowable values.
This function should be called on a pixel by pixel basis.
SUBROUTINE RELTAZ(YLAT,XLON,XNU,SZA,SUNZ,SUNA,SATZ,SATA,RELA,
$YLATS,XLONS,GHR,JDAY,DT,SDEC)
C
C**** GIVEN PIXEL LAT/LON AND SCAN ANGLE, COMPUTES SUB-SATELLLITE
C GROUND POINT FOR AVHRR.
C
C ** INPUT VARIABLES
C
C YLAT = THE LATITUDE EXPRESSED AS A VALUE BETWEEN -90 AND 90
C XLON = THE LONGITUDE EXPRESSED AS A VALUE BETWEEN 0 AND 360
C XNU = THE NEGATIVE OF THE SCAN ANGLE
C JDAY = THE JULIAN DAY, I.E. 001 IS JANUARY 1ST
C DT = THE EQUATORIAL CROSSING TIME, A FUNCTION OF SATELLITE
C SZA = THE ACTUAL SOLAR ZENITH ANGLE
C
C ** OUTPUT VARIABLES
C
C SUNZ = THE CALCULATED SOLAR ZENITH ANGLE
C SUNA = THE SOLAR AZIMUTH ANGLE
C SATZ = THE SATELLITE ZENITH ANGLE
C SATA = THE SATELLITE AZIMUTH ANGLE
C RELA = THE RELATIVE AZIMUTH ANGLE
C YLATS = THE SUB-SATELLITE LATITUDE
C XLONS = THE SUB-SATELLITE LONGITUDE
C GHR = THE GREENWICH MEAN TIME
C SDEC = THE SOLAR DECLINATION ANGLE, A FUNCTION OF JDAY.
C
INTEGER*4 JDAY
DOUBLE PRECISION RAD
DATA RAD /57.29578/
DATA RE /6378.0/
DATA HS /850.0/
C
C**** DATA CONSTANTS FOR SCENE
C
ROLL = 0.0
YAW = 0.0
TILT = 0.0
PERIOD = 101.94*60.0
EQT = DT
ICHECK=0
IC = 0
ICNT = 7
FEQT = 0.33
ANOMM = 50.
C
C**** OTHER IMPORTANT DATA
C
ALTCOR=(RE+HS)/RE
XINCLP=9.17
C
C**** PIXEL DATA
C
C**** THIS SUBROUTINE CALCULATES THE SOLAR DECLINATION.
C
CALL SUNAG1(RAD,JDAY,SDEC,TC)
C
C**** THIS SUBROUTINE CALCULATES THE SATELLITE ZENITH AND AZIMUTH
C ANGLES.
C
CALL SUBSAT(RAD,ALTCOR,XINCLP,PITCH,ROLL,YAW,TILT,XNU,
$XLON,YLAT,XLONS,YLATS,SATCA,SCZ)
SATA = SATCA
SATZ = SCZ
C
C**** THIS SUBROUTINE CALCULATES THE TIME ELAPSED SINCE EQUATORIAL
C CROSSING.
C
CALL SATEQT(RAD,PERIOD,XINCLP,XLONS,YLATS,
* XLON0S,TM)
99 GHR = (EQT - XLON0S/15.0) + TM/60.0
RLHR = GHR + XLON/15.0
C
C**** THIS SUBROUTINE CALCULATES THE SUN ZENITH AND AZIMUTH ANGLES.
C
CALL SUNAG2(RAD,GHR,SDEC,TC,XLON,YLAT,SUNZ,SUNA)
IF(SUNZ.GE.SZA-1..AND.SUNZ.LE.SZA+1.) THEN
RELA = ABS(SUNA-SATA)
C IF(RELA.GT.180) RELA = 360. - RELA
GO TO 100
ENDIF
IF(SUNZ.NE.SZA) THEN
IC = IC + 1
C
C**** MATCH UP CALCULATED SOLAR ZENITH WITH SATELLITE-GIVEN SOLAR
C ZENITH TO WITHIN +/- 1 DEGREE.
C
RELA = ABS(SUNA-SATA)
IF(ABS(SUNZ-SZA).LT.ANOMM) THEN
SAVSZ = SUNZ
SAVSA = SUNA
SAVRH = RLHR
SAVGH = GHR
SAVRE = RELA
ANOMM = ABS(SUNZ-SZA)
ENDIF
ENDIF
IF(SUNZ.GT.SZA) THEN
IF(IC.EQ.1) THEN
SUNZ = SZA
GO TO 100
ENDIF
SUNZ = SAVSZ
SUNA = SAVSA
RLHR = SAVRH
GHR = SAVGH
RELA = SAVRE
GO TO 100
ENDIF
FEQT = FEQT + 0.1
IHOUR = EQT
HOUR = IHOUR
IF(FEQT.LE.0.59) THEN
EQT = HOUR + FEQT
GO TO 99
ELSE
FEQT = FEQT - 0.6
EQT = HOUR + 1.0 + FEQT
GO TO 99
ENDIF
100 CONTINUE
RETURN
END
SUBROUTINE SATEQT(RAD,PERIOD,XINCLP,XLONS,YLATS,XLON0S,TM)
C
C GIVEN SATELLITE GROUND POINT (NADIR) LATITUDE AND LONGITUDE,
C COMPUTES TIME TRAVELLED (IN MINUTES) FROM LAST EQUATOR CROSSING
C AND LONGITUDE OF EQUATOR CROSSING.
C ONLY GOOD FOR ASCENDING NODE -- DO NOT EXCEED N OR S POLE.
C
DOUBLE PRECISION RAD
C
RINCLP = XINCLP/RAD
RLATS = YLATS/RAD
C
C CALCULATE degrees TRAVELED AND TRAVEL TIME
X = SIN(RLATS)/COS(RINCLP)
IF (X.GT.1.0) X = 1.0
IF (X.LT.-1.0) X = -1.0
RT = ASIN(X)
RT = ABS(RT)
DT = RT*RAD
TS = DT*PERIOD/360.0
C
C CALCULATE DELTA LONGITUDE
IF(YLATS.GE.0.) THEN
X = (SIN(RT)*SIN(RINCLP))/(COS(RLATS))
IF (X.GT.1.0) X = 1.0
IF (X.LT.-1.0) X = -1.0
RLAM=ASIN(X)
ELSE
X = (SIN(RT)*SIN(RINCLP))/(COS(RLATS))
IF (X.GT.1.0) X = 1.0
IF (X.LT.-1.0) X = -1.0
RLAM=-ASIN(X)
ENDIF
C
C CALCULATE EQUATOR CROSSING LONGITUDE
DLAM = RLAM*RAD
TTS = TS*4.167E-3
IF (YLATS .GE. 0.0)THEN
XLON0S = XLONS + DLAM + TTS
ELSE
XLON0S = XLONS - DLAM + TTS
ENDIF
C
C QUADRANT FIXUP
IF (XLON0S .GT. 180.0 .OR. XLON0S .LT. -180.0)THEN
IF (YLATS .GT. 0.0)THEN
XLON0S = XLON0S - 360.0
ELSE
XLON0S = XLON0S + 360.0
ENDIF
ENDIF
TM = TS/60.0
C
RETURN
END
SUBROUTINE SUBSAT(RAD,ALTCOR,XINCLP,PITCH,ROLL,YAW,TILT,
* XNU,XLONP,YLATP,XLONS,YLATS,SCA,SCZ)
C
C PURPOSE: THIS ROUTINE CALCULATES THE SATELLITE GROUND POINT
C LATITUDE AND LONGITUDE PLUS THE SPACECRAFT ZENITH
C AND THE TRUE BEARING FROM THE SATELLITE TO THE
C PIXEL.
C
C FORTRAN CALLING PROCEDURE:
C
C CALL SUBSAT (RAD,ALTCOR,XINCLP,PITCH,ROLL,YAW,TILT,XNU,
C XLONP,YLATP,XLONS,YLATS,SCA,SCZ)
C
C INPUT PARAMETERS:
C RAD = degrees/RADIANS CONVERSION FACTOR
C ALTCOR = ALTITUDE CORRECTION FACTOR= (R+H)/R = 1.149
C XINCLP = EQUATORIAL TRACKLINE INCLINATION = 9.28 DEGREES
C PITCH = PITCH ANGLE OF SPACECRAFT
C ROLL = ROLL ANGLE OF SPACECRAFT
C YAW = YAW OF SPACECRAFT
C TILT = SCAN MIRROR TILT ANGLE
C XNU = SCAN ANGLE TO PIXEL FOR ZERO TILT
C XLONP = LONGITUDE OF THE PIXEL
C YLATP = LATITUDE OF THE PIXEL
C
C OUTPUT PARAMETERS:
C YLATS = LATITUDE OF SPACECRAFT (SUB-SATELLITE POINT)
C XLONS = LONGITUDE OF SPACECRAFT
C SCA = AZIMUTH FROM NORTH SPACECRAFT TO PIXEL
C SCZ = SPACECRAFT ZENITH ANGLE
C
C CALLED BY: ATMO
C
C SUBROUTINES CALLED: CZSUBA
C
C VARIABLE LIST:
C RPSI = TRACKLINE RELATIVE BEARING TO PIXEL IN RADIANS
C RTHS = EFFECTIVE SCAN ANGLE MODIFIED BY TILT IN RADIANS
C RTH = EARTH SURFACE DISTANCE IN RADIANS FROM THE
C SATELLITE G. P. TO THE PIXEL.
C
C******************************************************************
C
C THE FIRST CALL IS TO CZSUBA FOR VARIOUS NAVIGATION PARAMETERS
DOUBLE PRECISION RAD
DATA RAD /57.29578/
CALL CZSUBA(RAD,ALTCOR,PITCH,ROLL,YAW,TILT,XNU,RPSI,RTHS,RTH)
RLATP = YLATP/RAD
RLONP = XLONP/RAD
RINCLP = XINCLP/RAD
C
C SET UP QUADRATIC VARIABLES FROM SCRIPPS ALGORITHM REPORT
C WILSON, ET AL. CZCS GEOLOCATION ALGORITHMS.
A = SIN(RLATP)-SIN(RTH)*SIN(RPSI)*SIN(RINCLP)
B = COS(RTH)
C = SIN(RTH) * COS(RPSI)
D = COS(RINCLP)**2
AB = A*B
XD = B*B+C*C
X = AB*AB-XD*(A*A-D*C*C)
XP = 0.0
IF (X.GT.1.0E-10)XP=SQRT(X)
C
C GET THE +/- VALUES FROM THE QUADRATIC AND
C IF TILT IS POSITIVE USE THE NEGATIVE QUADRATIC
C AND VICE VERSA
IF (TILT .LT. 0.0)THEN
XVAL = (AB+XP)/XD
ELSE
XVAL = (AB-XP)/XD
ENDIF
IF (XVAL.GT.1.0) XVAL = 1.0
IF (XVAL.LT.-1.0) XVAL = -1.0
RLATS = ASIN(XVAL)
YLATS = RLATS*RAD
C
C CALCULATE TRUE BEARING FROM NORTH TO PIXEL
C
X = SIN(RINCLP)/COS(RLATS)
IF (X.GT.1.0) X = 1.0
IF (X.LT.-1.0) X = -1.0
RINCLO = ASIN(X)
RPSIP = RPSI-RINCLO
TBEAR = RPSIP*RAD
IF (RPSI .LE. RINCLO)TBEAR = TBEAR+360.0
C
C RE-ORIENT SPACECRAFT AZIMUTH WITH RESPECT PIXEL
PI = 180.0/RAD
IF (RTH .NE. 0.0)THEN
COSCA = (SIN(RLATS)-COS(RTH)*SIN(RLATP))/(SIN(RTH)*COS(RLATP))
EPS = ABS(COSCA)
IF (EPS .GT. 1.0)THEN
IF (COSCA .LT. -1.0)COSCA=-1.0
IF (COSCA .GT. 1.0)COSCA=1.0
ENDIF
SCA = ACOS(COSCA)*RAD
IF (TBEAR .LE. 180.0)SCA = 360.0 - SCA
ELSE
IF (RPSI .LT. PI)SCA = RPSIP*RAD+180.0
IF (RPSI .GE. PI)SCA = RPSIP*RAD-180.0
ENDIF
C
C CALCULATE DELTA LONGITUDE
X = (COS(RTH)-SIN(RLATS)*SIN(RLATP))/(COS(RLATS)*COS(RLATP))
EPS = ABS(X)
IF (EPS .GE. 1.0)THEN
IF (X .LT. 0.0)X=-1.0
IF (X .GT. 1.0)X=1.0
ENDIF
DLAM = ACOS(X)
XDLAM = DLAM*RAD
C
C PIXEL LONGITUDE
IF (TBEAR .LE. 180.0)THEN
XVAL = XLONP - XDLAM
ELSE
XVAL = XLONP + XDLAM
ENDIF
IF (XVAL .LT. -180.0)THEN
XLONS = 360.0 + XVAL
ELSE IF (XVAL .GT. 180.0)THEN
XLONS = XVAL - 360.0
ELSE
XLONS = XVAL
ENDIF
C CALCULATE SPACECRAFT ZENITH ANGLE
R=ALTCOR*SIN(ABS(RTHS))
IF (R.GT.1.0) R = 1.0
IF (R.LT.-1.0) R = -1.0
TEMP=ASIN(R)
SCZ=TEMP*RAD
RETURN
END
C*****************************************************************
SUBROUTINE CZSUBA(RAD,ALTCOR,PITCH,ROLL,YAW,TILT,XNU,
* RPSI,RTHS,RTH)
C
C PURPOSE: THIS IS AN INTERNAL NAVIGATION SUBROUTINE CALLED
C MAINLY BY OTHER SUBROUTINES
C
C FORTRAN CALLING PROCEDURE:
C CALL CZSUBA(RAD,ALTCOR,PITCH,ROLL,YAW,TILT,XNU,RPSI,RTHS,RTH)
C
C INPUT PARAMETERS:
C RAD = DEGREES/RADIANS CONVERSION FACTOR
C ALTCOR = ALTITUDE CORRECTION FACTOR= (RE+HS)/RE=1.149
C ROLL,PITCH,YAW = ROLL,PITCH,YAW ANGLES OF SPACECRAFT
C TILT = MIRROR SHAFT TILT ANGLE
C XNU = SCAN ANGLE TO PIXEL FOR 0 TILT
C
C OUTPUT PARAMETERS:
C RPSI = TRACKLINE RELATIVE BEARING TO PIXEL IN RADIANS
C RTHS = EFFECTIVE SCAN ANGLE MODIFIED BY TILT
C RTH = EARTH SURFACE DISTANCE IN RADIANS G.P. TO PIXEL
C
C CALLED BY: SUBSAT
C
C **************************************************************
DOUBLE PRECISION RAD
DATA RAD /57.29578/
PID2 = ATAN(1.0)*2.0
PI = PID2*2.0
PI3D2 = PID2*3.0
RNU = XNU/RAD
T = TILT/2.0
RT = T/RAD
IF (TILT.NE.0.0) GOTO 1
RPSI = PID2
IF (XNU.LT.0.0) RPSI=PI3D2
RTHS = ABS(RNU)
GOTO 2
1 CONTINUE
C
2 CONTINUE
SINA = ALTCOR*SIN(RTHS)
IF(SINA.GT.1.0) SINA = 1.0
IF(SINA.LT.-1.0) SINA = -1.0
A = ASIN(SINA)
RTH = A-RTHS
RETURN
END
SUBROUTINE SUNAG1(RAD,IDAY,SDEC,TC)
C
C COMPUTES SOLAR DECLINATION AND TIME CORRECTION GIVEN
C JULIAN DAY. THIS IS TO BREAK UP THE SUNANG SUBROUTINE
C FOR BETTER COMPUTATIONAL EFFICIENCY. MUST BE COMBINED
C WITH SUNANG2 TO COMPUTE SOLAR ZENITH AND AZIMUTH.
C
C IDAY = JULIAN DAY
C SDEC = SOLAR DECLINATION ANGLE IN DEGREES
C THEZ = THETA ZERO ORBITAL POSITION IN DEGREES
C TC = TIME CORRECTION
C
C COMPUTE SOLAR DECLINATION ANGLE
DOUBLE PRECISION RAD
THEZ = 360.0*(IDAY-1)/365.0
RTHEZ = THEZ/RAD
SDEC = 0.396372-22.91327*COS(RTHEZ) + 4.02543*SIN(RTHEZ)
* - 0.387205*COS(2.0*RTHEZ) + 0.051967*SIN(2.0*RTHEZ)
* - 0.154527*COS(3.0*RTHEZ) + 0.084798*SIN(3.0*RTHEZ)
RSDEC = SDEC/RAD
C
C TIME CORRECTION FOR SOLAR HOUR ANGLE, AND SOLAR HOUR ANGLE
TC = 0.004297 + 0.107029*COS(RTHEZ) - 1.837877*SIN(RTHEZ)
* - 0.837378*COS(2.0*RTHEZ) - 2.340475*SIN(2.0*RTHEZ)
C
RETURN
END
C ***************************************************************
SUBROUTINE SUNAG2(RAD,HR,SDEC,TC,XLON,YLAT,SUNZ,SUNA)
C
C COMPUTES SUN AZIMUTH AND ZENITH ANGLES FOR A GIVEN
C TIME, DATE, LATITUDE AND LONGITUDE. THIS PROGRAM
C IS FROM THE NMFS ELAS COMPUTER CODE. MODIFIED FOR
C STANDARD COORDINATES (W LONG. NEGATIVE), AND TO CORRECT
C FOR DATELINE PROBLEM. THIS IS PART 2 OF THE FORMER
C SUNANG SUBROUTINE. THIS REQUIRES DECLINATION AND TIME
C CORRECTION COMPUTED FROM SUNANG1.
C
C TIME = TIME OF DAY IN SECONDS
C YLAT = LATITUDE OF PIXEL
C XLON = LONGITUDE OF PIXEL
C SUNZ = SOLAR ZENITH IN DEGREES
C SUNA = SOLAR AZIMUTH ANGLE (FROM NORTH) IN DEGREES
C SDEC = SOLAR DECLINATION ANGLE IN DEGREES
C THEZ = THETA ZERO ORBITAL POSITION IN DEGREES
C TC = TIME CORRECTION
C XHA = SOLAR HOUR ANGLE IN DEGREES
C
DOUBLE PRECISION RAD
RSDEC = SDEC/RAD
C
C SOLAR HOUR ANGLE
XHA = (HR-12.0)*15.0 + XLON + TC
IF (XHA .GT. 180.0)XHA = XHA-360.0
IF (XHA .LT. -180.0)XHA = XHA+360.0
RLAT = YLAT/RAD
RLON = XLON/RAD
RHA = XHA/RAD
C
C SUN ZENITH
COSTMP = SIN(RLAT)*SIN(RSDEC) +
* COS(RLAT)*COS(RSDEC)*COS(RHA)
EPS = ABS(COSTMP)
IF (EPS .GT. 1.1)THEN
ELSE IF (EPS .GT. 1.0)THEN
IF (COSTMP .GT. 0.0)COSTMP=1.0
IF (COSTMP .LT. 0.0)COSTMP=-1.0
ENDIF
RSUNZ = ACOS(COSTMP)
C
C SUN AZIMUTH
COSAZ = (SIN(RSDEC)-SIN(RLAT)*COS(RSUNZ))/(COS(RLAT)*SIN(RSUNZ))
EPS = ABS(COSAZ)
IF (EPS .GT. 1.1)THEN
ELSE IF (EPS .GT. 1.0)THEN
IF (COSAZ .GT. 0.0)COSAZ=1.0
IF (COSAZ .LT. 0.0)COSAZ=-1.0
ENDIF
RSUNA = ACOS(COSAZ)
IF (XHA .GT. 0.0)RSUNA=360.0/RAD-RSUNA
C
C CONVERT TO DEGREES
SUNZ = RSUNZ*RAD
SUNA = RSUNA*RAD
C
RETURN
END
| Previous Section | Top of Page | Next Section |