PROGRAM FRESNO C C Compute the position of the sun for Fresno on the 15th of the C month from April to September. C REAL MIN,LAT,LONG,MONTH,MONTHS(6) INTEGER*2 DATE(6) C DATA DATE / 106,136,167,197,228,259 / DATA MONTHS / 'APR ','MAY ','JUN ','JUL ','AUG ','SEP ' / C IYEAR = 1984 YEAR = IYEAR IDAY = 228 DAY = IDAY MIN = 0. SEC = 0. DASVTM= 0. LAT = 36.767 LONG = 119.717 ZONE = 8. C DO 1000 N=1,6 IF (N .EQ. 1) THEN DASVTM = 0. ELSE DASVTM = 1. ENDIF C IDAY = DATE(N) DAY = IDAY MONTH = MONTHS(N) WRITE (6,10) IYEAR,IDAY,MONTH,LAT,LONG 10 FORMAT (10X'P O S I T I O N O F T H E S U N'// 1 10X'YEAR:'I4/ 2 10X'DAY: 'I3' ('A4'15)'/ 3 10X'LATITUDE = 'F7.3' DEGREES NORTH'/ 4 10X'LONGITUDE = 'F7.3' DEGREES WEST OF GREENWICH'//) IF (DASVTM .NE. 0.) WRITE (6,20) 20 FORMAT (10X'Daylight Savings Time in affect.'//) WRITE (6,30) 30 FORMAT (10X'HOUR'10X'AZIMUTH'10X'ELEVATION'//) C DO 100 I=0,24 HR=I CALL SUNAE(YEAR,DAY,HR,MIN,SEC,ZONE,DASVTM,LAT,LONG,A,E) IZ='AM' J=I IF(HR.GT.12.) J=J-12 IF(HR.GE.12.) IZ='PM' IF (E .GE. -10.) WRITE(6,40) J,IZ,A,E 40 FORMAT(9X,I2,1X,A2,F16.2,F18.2) 100 CONTINUE C WRITE(6,110) 110 FORMAT(1H1) C 1000 CONTINUE C CALL EXIT END C************************************************************************* SUBROUTINE SUNAE(YEAR,DAY,HR,MIN,SEC,ZONE,DASVTM,LAT,LONG,SA,SE) C C THIS SUBROUTINE CALCULATES THE LOCAL AZIMUTH AND ELEVATION OF C THE SUN AT A SPECIFIED LOCATION AND TIME USING AN APPROXIMATION C TO THE EQUATIONS USED TO GENERATE THE NAUTICAL ALMANAC. C C INPUT PARAMETERS C YEAR - THE YEAR NUMBER (E.G., 1977). C DAY - THE DAY NUMBER OF THE YEAR STARTING WITH 1 FOR C JANUARY 1, EXCEPT IN LEAP YEARS WHEN 1 SHOULD BE C SUBTRACTED FROM THE DAY NUMBER BEFORE MARCH 1. C HR,MIN,SEC - THE TIME OF DAY. C ZONE - THE LOCAL INTERNATIONAL ZONE TIME (E.G.: P.S.T. = 8). C DASVTM - =1 IF DAYLIGHT SAVINGS TIME IN EFFECT, ELSE = 0. C LAT - THE LOCAL LATITUDE IN DEGREES (NORTH IS POSITIVE). C LONG - THE LOCAL LONGITUDE IN DEGREES WEST OF GREENWICH. C C OUTPUT PARAMETERS C SA - AZIMUTHAL ANGLE OF THE SUN (POSITIVE IS EAST OF SOUTH) C SE - ELEVATION OF THE SUN. C C------------------------------------------------------------------------- C REAL MIN, LAT, LONG DATA TWOPI, RAD /6.2831853, 0.017453293/ C DELYR = YEAR - 1980. LEAP = IFIX(DELYR/4.) T = HR + (MIN + SEC/60.) / 60. + ZONE - DASVTM TIME = DELYR*365. + LEAP + DAY - 1. + T/24. IF (DELYR .EQ. LEAP*4.) TIME = TIME - 1. IF ((DELYR.LT.0.) .AND. (DELYR.NE.LEAP*4.)) TIME = TIME - 1. THETA = (360.*TIME/365.25) * RAD G = -0.031271 - 4.53963E-7*TIME + THETA EL = 4.900968+3.6747E-7*TIME+(0.033434-2.3E-9*TIME)*SIN(G) # +0.000349*SIN(2.*G) + THETA EPS = 0.409140 - 6.2149E-9*TIME SEL = SIN(EL) A1 = SEL * COS(EPS) A2 = COS(EL) RA = ATAN2(A1,A2) IF (RA .LT. 0.) RA = RA + TWOPI DECL = ASIN(SEL*SIN(EPS)) ST = 1.759335 + TWOPI*(TIME/365.25-DELYR) + 3.694E-7*TIME IF (ST .GE. TWOPI) ST = ST - TWOPI S = ST + (T*15. - LONG) * RAD IF (S .GE. TWOPI) S = S - TWOPI H = RA - S PHI = LAT * RAD TEMP = ASIN(SIN(PHI)*SIN(DECL)+COS(PHI)*COS(DECL)*COS(H)) SE = TEMP SA = ASIN(COS(DECL)*SIN(H)/COS(SE))/RAD IF (SIN(SE) .GE. SIN(DECL)/SIN(PHI)) GO TO 10 IF (SA .LT. 0.) SA = SA + 360. SA = 180. - SA 10 SE = SE / RAD RETURN END C*********************************************************************** FUNCTION ASIN(XX) X = XX IF (XX.LT.0.) X=-X Y=3.14159265/2.-SQRT(1.-X)*(1.5707963050+X*(-.2145988016 : +X*(.0889789874+X*(-.0501743046+X*(.0308918810+X*( : -.0170881256+X*(.0066700901+X*(-.0012624911)))))))) IF(XX.LT.0.)Y=-Y ASIN=Y RETURN END