parameter (mxtim=100) real tary(mxtim) real*8 cozena, dlat, dlon, dday, dtime, val REAL lon, lat, musun real jd, month, day, year, iday, hour, minute, sec Write(*,*) '--------------------------------------------------' Write(*,*) ' PROGRAM SUN_Position ver. 1.0 ' Write(*,*) ' Oblicza polozenie Slonca dla czasu uniwersalnego ' Write(*,*) ' Autor Krzysztof Markowicz ' Write(*,*) ' e-mail: kmark@igf.fuw.edu.pl ' Write(*,*) ' www.igf.fuw.edu.pl/meteo/stacja/ ' Write(*,*) '--------------------------------------------------' Write(*,*) ' Podaj szerokosc geograficzna:' Read(*,*) lat Write(*,*) ' Podaj dlugosc geograficzna:' Read(*,*) lon Write(*,*) ' Podaj rok: ' Read(*,*) year Write(*,*) ' Podaj miesiac:' Read(*,*) month Write(*,*) ' Podaj dzien:' Read(*,*) day CALL date2jd(year,month,day,jd) day=jd Write(*,*) ' Podaj pelna godzine w czasie uniwersalnym:' Read(*,*) hour Write(*,*) ' Podaj ilosc minut:' Read(*,*) minute Write(*,*) ' Podaj ilosc sekund:' Read(*,*) sec time=hour+minute/60+sec/3600 pi=4.*atan(1.) dtor=pi/180. xnull=-999. ntim=24*10 5 continue c print*, ' year day lat lon ', year, day, lat, lon, time c call sunae(year,day,time,lat,lon,a,e) call sunae1(year,day,time,lat,lon,az,el) c CALL pro(day,time,-lon,lat,zenith,azimuth) c call SUNANG( DAY, TIME, LON, LAT, MUSUN, RSUN ) print*,'**********************************************' print*,'Polozenie slonca:' print*,'Kat elewacyjny oraz azymut Slonca' 105 format (f6.1, f6.1) write(6,105) el, az dday=day dtime=time dlat=lat dlon=lon val=acos(COZENA(dDAY,dtime,DLAT,DLON))/dtor c write(6,56) time, acos(musun)/dtor,zenith,val 56 format(1x,4f12.6) 10 continue STOP END subroutine date2jd(year,month,day,jd) real year, month, day, jd, Feb, x x=mod(year,4.0) if (x.eq.0) then Feb=29 else Feb=28 end if if (month.eq.1) then jd=day end if if (month.eq.2) then jd=31+day end if if (month.eq.3) then jd=31+Feb+day end if if (month.eq.4) then jd=31+Feb+31+day end if if (month.eq.5) then jd=31+Feb+31+30+day end if if (month.eq.6) then jd=31+Feb+31+30+31+day end if if (month.eq.7) then jd=31+Feb+31+30+31+30+day end if if (month.eq.8) then jd=31+Feb+31+30+31+30+31+day end if if (month.eq.9) then jd=31+Feb+31+30+31+30+31+31+day end if if (month.eq.10) then jd=31+Feb+31+30+31+30+31+31+30+day end if if (month.eq.11) then jd=31+Feb+31+30+31+30+31+31+30+31+day end if if (month.eq.12) then jd=31+Feb+31+30+31+30+31+31+30+31+30+day end if return end subroutine none C W. H. Wilson, "Solar ephemeris algorithm," Reference 80-13, 70 pp., C Scripps Institution C of Oceanography, University of California at San Diego, C La Jolla, California, 1980. C http://ice.cor.epa.gov/projects/ipw/man/pages/sunang.html return end SUBROUTINE pro(day,time,lon,lat,zenith,azimuth) REAL nday(74), eqt(74), dec(74) REAL lon, lat REAL latsun, lonsun DATA nday/1.0, 6.0, 11.0, 16.0, 21.0, 26.0, 31.0, 36.0, 41.0, $ 46.0, 51.0, 56.0, 61.0, 66.0, 71.0, 76.0, 81.0, 86.0, 91.0, $ 96.0, 101.0, 106.0, 111.0, 116.0, 121.0, 126.0, 131.0, 136.0, $ 141.0, 146.0, 151.0, 156.0, 161.0, 166.0, 171.0, 176.0, $ 181.0, 186.0, 191.0, 196.0, 201.0, 206.0, 211.0, 216.0, $ 221.0, 226.0, 231.0, 236.0, 241.0, 246.0, 251.0, 256.0, $ 261.0, 266.0, 271.0, 276.0, 281.0, 286.0, 291.0, 296.0, $ 301.0, 306.0, 311.0, 316.0, 321.0, 326.0, 331.0, 336.0, $ 341.0, 346.0, 351.0, 356.0, 361.0, 366.0/ DATA eqt/-3.23, -5.49, -7.60, -9.48, -11.09, -12.39, -13.34, $ -13.95, -14.23, -14.19, -13.85, -13.22, -12.35, -11.26, $ -10.01, -8.64, -7.18, -5.67, -4.16, -2.69, -1.29, -0.02, $ 1.10, 2.05, 2.80, 3.33, 3.63, 3.68, 3.49, 3.09, 2.48, 1.71, $ 0.79, -0.24, -1.33, -2.41, -3.45, -4.39, -5.20, -5.84, -6.28, $ -6.49, -6.44, -6.15, -5.60, -4.82, -3.81, -2.60, -1.19, 0.36, $ 2.03, 3.76, 5.54, 7.31, 9.04, 10.69, 12.20, 13.53, 14.65, $ 15.52, 16.12, 16.41, 16.36, 15.95, 15.19, 14.09, 12.67, $ 10.93, 8.93, 6.70, 4.32, 1.86, -0.62, -3.23/ DATA dec/-23.06, -22.57, -21.91, -21.06, -20.05, -18.88, -17.57, $ -16.13, -14.57, -12.91, -11.16, -9.34, -7.46, -5.54, -3.59, $ -1.62, 0.36, 2.33, 4.28, 6.19, 8.06, 9.88, 11.62, 13.29, $ 14.87, 16.34, 17.70, 18.94, 20.04, 21.00, 21.81, 22.47, $ 22.95, 23.28, 23.43, 23.40, 23.21, 22.85, 22.32, 21.63, $ 20.79, 19.80, 18.67, 17.42, 16.05, 14.57, 13.00, 11.33, 9.60, $ 7.80, 5.95, 4.06, 2.13, 0.19, -1.75, -3.69, -5.62, -7.51, $ -9.36, -11.16, -12.88, -14.53, -16.07, -17.50, -18.81, $ -19.98, -20.99, -21.85, -22.52, -23.02, -23.33, -23.44, $ -23.35, -23.06/ C INPUT: C day Julian day C time Universal Time in hours C lat geographic latitude of point on earth's surface (degrees) C lon geographic longitude of point on earth's surface (degrees) C OUTPUT: C zenith solar zenith angle (degrees) C azimuth solar azimuth (degrees) C Azimuth is measured counter-clockwise from due south C PROCEDURE: C 1. Calculate the subsolar point latitude and longitude, based on C DAY and TIME. Since each year is 365.25 days long the exact C value of the declination angle changes from year to year. For C precise values consult THE AMERICAN EPHEMERIS AND NAUTICAL C ALMANAC published yearly by the U.S. govt. printing office. C The subsolar coordinates used in this code were provided by C a program written by Jeff Dozier. C C 2. Given the subsolar latitude and longitude, spherical geometry is C used to find the solar zenith, azimuth and flux multiplier. C C AUTHOR: Paul Ricchiazzi 23oct92 C Earth Space Research Group, UCSB C C SOURCE: This routine was adapted from the LOWTRAN (vers7) subroutine C SUBSOL.F. Some corrections were required C C Translated back to FORTRAN by P. J. Flatau pi = 4.*atan(1.) dtor = pi/180. dd = mod(real(day)+time/24.-1.,365.25) + 1. DO 10 i = 2, 74 IF (dd .LE. nday(i)) THEN k = i - 1 eqtime = eqt(k) + (eqt(k+1)-eqt(k))* (dd-nday(k))/ $ (nday(k+1)-nday(k)) eqtime = eqtime/60. decang = dec(k) + (dec(k+1)-dec(k))* (dd-nday(k))/ $ (nday(k+1)-nday(k)) GO TO 20 END IF 10 CONTINUE 20 CONTINUE latsun = decang lonsun = -15.* (time-12.+eqtime) t0 = (90.-lat)*dtor t1 = (90.-latsun)*dtor p0 = lon*dtor p1 = lonsun*dtor zz = cos(t0)*cos(t1) + sin(t0)*sin(t1)*cos(p1-p0) yy = sin(t1)*sin(p1-p0) xx = cos(t0)*sin(t1)*cos(p1-p0) - sin(t0)*cos(t1) azimuth = -atan2(yy,xx)/dtor zenith = acos(zz)/dtor rsun = 1. - 0.01673*cos(.9856* (dd-2.)*dtor) solfac = zz/rsun**2 RETURN END SUBROUTINE SUNANG( DAY, TIME, LONG, LAT, MUSUN, RSUN ) C from Wiscombe (ATRAD) C C CALCULATE SOLAR ZENITH ANGLE AND EARTH-SUN DISTANCE FROM C GEOGRAPHICAL AND TEMPORAL INFORMATION C C INPUT : DAY = DAY OF YEAR ( JAN 1 IS DAY 1 ) C TIME = GREENWICH TIME ( HOURS ) C LAT = LATITUDE ( DEGREES ) C LONG = LONGITUDE ( DEGREES ) C C OUTPUT : MUSUN = COSINE OF SOLAR ZENITH ANGLE C RSUN = SQUARE OF RATIO OF MEAN EARTH-SUN DISTANCE C TO EARTH-SUN DISTANCE AT *DAY* C C *** ARGUMENTS REAL DAY, TIME, LONG, LAT, MUSUN, RSUN C C *** LOCAL VARIABLES C C COSDLT COSINE OF DECLINATION ANGLE C D ANGLE MEASURED FROM PERIHELION (WHICH IS C APPROXIMATED AS MIDNITE JAN 1), IN RADIANS C DPY DAYS PER YEAR C EPSILN ECCENTRICITY OF EARTH ORBIT C H HOUR ANGLE (DEGREES) C M MERIDIAN PASSAGE (HOURS) C PI 3.14159... C SINDLT SINE OF DECLINATION ANGLE C SINOB SINE OF OBLIQUITY OF ECLIPTIC C REAL COSDLT, D, DPY, EPSILN, H, M, L, PI, SIGMA, SINDLT, SINOB DATA EPSILN / .016733 /, SINOB / 0.3978 /, DPY / 365.242 / C C PI = 2. * ASIN( 1.0 ) D = 2. * PI * ( DAY - 1. ) / DPY M = 12. + 0.12357 * SIN( D ) - 0.004289 * COS( D ) $ + 0.153809 * SIN( 2. * D ) + 0.060783 * COS( 2. * D ) H = 15. * ( TIME - M ) - LONG L = 279.9348 * (PI/180.) + D SIGMA = L + ( PI/180. ) * ( 0.4087 * SIN( L ) + 1.8724 * COS( L ) $ - 0.0182 * SIN( 2. * L ) + 0.0083 * COS( 2. * L ) ) SINDLT = SINOB * SIN( SIGMA ) COSDLT = SQRT( 1. - SINDLT**2 ) MUSUN = SINDLT * SIN( PI/180. * LAT ) $ + COSDLT * COS( PI/180. * LAT ) * COS( ( PI/180. ) * H ) RSUN = 1. / ( 1. - EPSILN * COS( D ) )**2 C RETURN END subroutine sunae1(year,day,hour,lat,long,az,el) implicit real(a-z) C Purpose: C Calculates azimuth and elevation of sun C C References: C (a) Michalsky, J. J., 1988, The Astronomical Almanac's algorithm for C approximate solar position (1950-2050), Solar Energy, 227---235, 1988 C C (b) Spencer, J. W., 1989, Comments on The Astronomical C Almanac's algorithm for approximate solar position (1950-2050) C Solar Energy, 42, 353 C C Input: C year - the year number (e.g. 1977) C day - the day number of the year starting with 1 for C January 1 C time - decimal time. E.g. 22.89 (8.30am eastern daylight time is C equal to 8.5+5(hours west of Greenwich) -1 (for daylight savings C time correction C lat - local latitude in degrees (north is positive) C lon - local longitude (east of Greenwich is positive. C i.e. Honolulu is 15.3, -157.8) C C Output: C a - azimuth angle of the sun (measured east from north 0 to 360) C e - elevation of the sun C C Spencer correction introduced and 3 lines of Michalsky code C commented out (after calculation of az) C pi=4.*atan(1.) twopi=2.*pi rad=pi/180. C get the current Julian date delta=year-1949. leap=aint(delta/4.) jd=32916.5+delta*365.+leap+day+hour/24. C calculate ecliptic coordinates time=jd-51545.0 C force mean longitude between 0 and 360 degs mnlong=280.460+0.9856474*time mnlong=mod(mnlong,360.) if(mnlong.lt.0) mnlong=mnlong+360. C mean anomaly in radians between 0, 2*pi mnanom=357.528+0.9856003*time mnanom=mod(mnanom,360.) if(mnanom.lt.0.) mnanom=mnanom+360. mnanom=mnanom*rad C compute ecliptic longitude and obliquity of ecliptic C eclong=mnlong+1.915*(mnanom)+0.20*sin(2.*mnanom) eclong=mnlong+1.915*sin(mnanom)+0.020*sin(2.*mnanom) eclong=mod(eclong, 360.) if(eclong.lt.0) eclong=eclong+360. oblqec=23.429-0.0000004*time eclong=eclong*rad oblqec=oblqec*rad C calculate right ascention and declination num=cos(oblqec)*sin(eclong) den=cos(eclong) ra=atan(num/den) C force ra between 0 and 2*pi if(den.lt.0.) then ra=ra+pi elseif(num.lt.0.) then ra=ra+twopi endif C dec in radians dec=asin(sin(oblqec)*sin(eclong)) C calculate Greenwich mean sidereal time in hours gmst=6.697375+0.0657098242*time+hour C hour not changed to sidereal sine "time" includes the fractional day gmst=mod(gmst,24.) if(gmst.lt.0.) gmst=gmst+24. C calculate local mean sidereal time in radians lmst=gmst+long/15. lmst=mod(lmst,24.) if(lmst.lt.0.) lmst=lmst+24. lmst=lmst*15.*rad C calculate hour angle in radians between -pi, pi ha = lmst -ra if(ha.lt.-pi) ha=ha+twopi if(ha.gt.pi) ha=ha-twopi lat=lat*rad C calculate azimuth and elevation el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha)) az=asin(-cos(dec)*sin(ha)/cos(el)) C add J. W. Spencer code (next 5 lines) if(sin(dec) - sin(el)*sin(lat) .ge. 0. ) then if(sin(az) .lt. 0.) az=az+twopi else az=pi-az endif C end Spencer's corrections ccc elc=asin(sin(dec)/sin(lat)) ccc if(el.ge.elc) az=pi-az ccc if(el.le.elc .and. ha.gt.0.) az=twopi+az C this puts azimuth between 0 and 2*pi radians C calculate refraction correction for US stand. atm. el=el/rad if(el.gt.-0.56) then refrac=3.51561*(0.1594+0.0196*el+0.00002*el**2)/ $ (1.+0.505*el+0.0845*el**2) else refrac=0.56 endif el=el+refrac az=az/rad lat=lat/rad C return end subroutine sunae(year,day,t,lat,long,a,e) real lat, long C Purpose: C Calculates azimuth and elevation of sun C C References: C (a) Walraven, R., 1978, Calculating the position of the sun, C Solar Energy, 20, 393---397 C (b) Walraven, R., 1979, Erratum: Calculating the position of the sun, C Solar Energy, 22, 195 C Input: 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 subtracted C from the day number before March 1 C time - decimal time. E.g. 22.89 C lat - local latitude in degrees (north is positive) C lon - local longitude in degrees west of Greenwich e.g. 121.758 C Honolulu is (15.3, 157.8) C Output: C a - azimuth angle of the sun (postive is east of south) between 0-180 C e - elevation of the sun C pi=4.*atan(1.) twopi=2.*pi rad=pi/180. delyr=year-1980. leap=ifix(delyr/4.) C corrected in Walraven, 1979; and by PJF to include decimal time C 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 C corrected Walraven, 1979 s=st+(t*15.-long)*rad if(s.ge.twopi) s=s-twopi h=ra-s phi=lat*rad e=asin(sin(phi)*sin(decl)+cos(phi)*cos(decl)*cos(h)) a=asin(cos(decl)*sin(h)/cos(e))/rad if(sin(e) .ge. sin(decl)/sin(phi)) go to 10 if(a .lt. 0.) a=a+360. a=180.-a 10 e=e/rad return end subroutine sunang1(iday,hr,rad,xlon,ylat,sunz) C from Mobley HYDROLIGHT 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), to correct c for dateline problem, and to correct coefficients (taken c from Iqbal, 1983, An Introduction to Solar Radiation). c Watson Gregg, Research and Data Systems, Corp. c c iday = julian day c time = time of day in seconds c ylat = latitude of pixel c xlon = longitude of pixel c sunz = solar zenith 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 c c Compute solar declination angle 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.342824*sin(2.0*rthez) 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 write(6,*)'Error in acos argument in sun zenith' write(6,*)costmp 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 Sun azimuth sintmp = sin(abs(rha))*cos(rsdec)/sin(rsunz) rsuna = asin(sintmp) if(ylat .gt. sdec) rsuna = 180.0 / rad - rsuna if(xha .gt. 0.) rsuna = 360.0/rad - rsuna c c Convert to degrees suna = rsuna * rad c write(6,*),"The sun's computed azimuth is",suna,' degrees' sunz = rsunz*rad c return end real*8 function ZENITH (YR, MONTH, DAY, HOURZ, LAT, LON) c-------------------------------------------------------------------------- C Calculates solar zenith angle, returning the cosine. c The routine then calls the function c cozena which returns the cosine of the zenith angle c c Input: c yr Year (e.g., 92) [integer] c month,day Month (e.g., 7), day of month (e.g., 15) [integer] c hourz Decimal hour (GMT), e.g., 13.2 [real*8] c tlat,tlon Latitude (positive in northern hemisphere), longitude c (0 to +- 180 OR 0 to 360, positive west of Prime Mer.) c [real*8] c c Functions called: c COZENA(DAY,HOUR,LAT,TLON) c-------------------------------------------------------------------------- implicit none INTEGER MONTH, YR, DAY REAL*8 LAT,lon,jday,hourz,cozena logical leap c-------------------------------------------------------------------------- c Leap year: 'leap'=.true. if a leap year. if (mod(yr+1900,4) .eq. 0) then leap = .true. else leap = .false. end if C FIND TIME IN JULIAN DAYS: C from STREAMER IF (MONTH.EQ.1) jday = 0. IF (MONTH.EQ.2) jday = 31. IF (MONTH.EQ.3) jday = 59. IF (MONTH.EQ.4) jday = 90. IF (MONTH.EQ.5) jday = 120. IF (MONTH.EQ.6) jday = 151. IF (MONTH.EQ.7) jday = 181. IF (MONTH.EQ.8) jday = 212. IF (MONTH.EQ.9) jday = 243. IF (MONTH.EQ.10) jday = 273. IF (MONTH.EQ.11) jday = 304. IF (MONTH.EQ.12) jday = 334. IF (MONTH .GT. 2 .AND. LEAP) jday = jday + 1. jday = jday + DAY zenith = cozena(jday,hourz,lat,lon) return end c***************************************************************************** DOUBLE PRECISION FUNCTION COZENA(DAY,HOUR,DLAT,DLON) C+---------------------------------------------------------------------+ C! ! C! REFERENCE: Woolf, H.M., NASA TM X-1646. ! C! ON THE COMPUTATION OF SOLAR ELEVATION ! C! ANGLES AND THE DETERMINATION OF SUNRISE ! C! AND SUNSET TIMES. ! C+---------------------------------------------------------------------+ C! I N P U T V A R I A B L E S: ! C+---------------------------------------------------------------------+ C! DAY DAYS OF YEAR (JAN. 1 IS DAY 1) ! C! HOUR HOURS OF GREENWICH TIME ! C! DLON DEGREE LONGITUDE OF GRID POINT, ! C! COUNTED POSITIVE WEST OF GREENWICH ! C! DLAT DEGREE LATITUDE OF GRID POINT, ! C! COUNTED POSITIVE NORTHERN HEMISPHERE ! C+---------------------------------------------------------------------+ C! O U T P U T V A R I A B L E S: ! C+---------------------------------------------------------------------+ C! COZENA COSINE OF ZENITH ANGLE ( Reference -- EQ. 1.1 ) ! C+---------------------------------------------------------------------+ C! D A T A V A R I A B L E S: ! C+---------------------------------------------------------------------+ C! EPSILN ECCENTRICITY OF EARTH ORBIT ! C! SINOB SINE OF OBLIQUITY OF ECLIPTIC ! C! DPY DAYS PER YEAR (365.242) ! C! DPH DEGREES PER HOUR (360./24) ! C+---------------------------------------------------------------------+ C! I N T E R N A L V A R I A B L E S: ! C+---------------------------------------------------------------------+ C! DPR DEGREE/RADIAN (*57.29578) ! C! RPD RADIAN/DEGREE (*.01745329) ! C! DANG ANGLE MEASURED FROM PERIHELION, IN RADIANS ! C! WHICH IS TAKEN AS MIDNIGHT JAN. 1 ! C! HOMP HOURS OF MERIDIAN PASSAGE OR TRUE SOLAR NOON ! C! ( Reference -- EQ. 1.6 ) ! C! HANG HOUR ANGLE, A MEASURE OF THE LONGITUDINAL DISTANCE ! C! TO THE SUN FROM THE POINT CALCULATED ( EQ. 1.5 ) ! C! SINDLT SINE OF DECLINATION ANGLE ( EQ. 1.2 ) ! C! COSDLT COSINE OF DECLINATION ANGLE ! C! SIGMA Reference -- EQ. 1.3A ! C! ANG Reference -- EQ. 1.3B ! C+---------------------------------------------------------------------+ IMPLICIT REAL*8 (A-H,O-Z) DATA EPSILN/.016733/,SINOB/.3978/,DPY/365.242/,DPH/15.0/ C from STREAMER PI = 3.1415927 RPD = PI/180.0 DPR = 1.0/RPD DANG = 2.0*PI* (DAY-1.0)/DPY HOMP = 12.0 + 0.123570*SIN(DANG) - 0.004289*COS(DANG) + & 0.153809*SIN(2.0*DANG) + 0.060783*COS(2.0*DANG) HANG = DPH* (HOUR-HOMP) - DLON ANG = 279.9348*RPD + DANG SIGMA = (ANG*DPR+0.4087*SIN(ANG)+1.8724*COS(ANG)- & 0.0182*SIN(2.0*ANG)+0.0083*COS(2.0*ANG))*RPD SINDLT = SINOB*SIN(SIGMA) COSDLT = SQRT(1.0-SINDLT**2) COZENA = SINDLT*SIN(RPD*DLAT) + COSDLT*COS(RPD*DLAT)*COS(RPD*HANG) RETURN END