PROGRAM solpositiondblfunc
C
C Is based on the program of Drs J. Michalski
C and L. Harrison, using the
C Astronomical Almanac algorithm
C and more specifically the "sunae" subroutine
C which is Dr W. Winscombe's routine incorporating
C modifications and fixes. The refraction adjustment
C is the NOAA method.
C
C The calculations are broken down into simpler
C steps to suit a simple f77 interpreter.
C It needs input checking, further testing etc and is
C working at least for Australian coordinates, 2000-2050.
C It was an exercise to test the interpreter.
C Errors are almost certainly all mine :)
C S. Oliver Sept 2015
C
C Input :
C
C latitude, longitude S = -ve, E = +ve
C Day,month,year eg 1,12,2015
C = 1 Dec 2015 (range 2000 - 2050)
C Time zone offset Hours ahead of GMT
C eg Sydney is 10.0.
C
C Output :
C
C Day Day of year
C
C Hr Hours from midnight
C
C El Elevation ( degrees up from horizon )
C Az Azimuth ( degrees cw from north )
C
C
DOUBLE az, el, hour, zone, lat, long, twil
DOUBLE dec, den, eclong, gmst, ha, jd, lmst
DOUBLE mnanom, mnlong, num, oblqec, pi, ra
DOUBLE rpd, refrac, time, twopi, temp, eltan
INTEGER*4 i, year, day, jtemp
INTEGER*4 delta, leap
INTEGER*2 dy, mnth, yr, rem
INTEGER*2 leapday, yearday
C
FORMAT (F4.1,X,F9.4,X,F9.4,/)
C
pi = 3.14159265359
twopi = 2.0 * pi
rpd = 0.01745329252
C
C
year = 2000
day = 1
hour = 0.0
temp = 360.0
lat = 0.0
long = 0.0
C
C
PRINT *, "Latitude and longitude ?"
READ lat, long
PRINT *, "Day, month, year ?"
READ dy, mnth, yr
PRINT *, "What is the time zone offset ?"
READ zone
C
C
yearday = dy + ( mnth - 1 ) * 31
C
C jan feb mar apr may jun jul aug sep oct nov dec
C 31 28 31 30 31 30 31 31 30 31 30 31
C
C feb, apr, jun, sep, nov
IF ( mnth > 2 ) yearday = yearday - 3
IF ( mnth > 4 ) yearday = yearday - 1
IF ( mnth > 6 ) yearday = yearday - 1
IF ( mnth > 9 ) yearday = yearday - 1
IF ( mnth > 11 ) yearday = yearday - 1
C
leapday = 0
rem = yr - ( yr / 4 ) * 4
IF ( rem = 0 ) leapday = 1
IF ( mnth > 2 ) yearday = yearday + leapday
C
PRINT *, "Day of year is ", yearday
PRINT
PRINT *, "Hr El Az"
day = yearday
year = yr
C
hour = 0.0 - zone
C
DO 10 i = 1, 24
hour = hour + 1.0
temp = 24.0
C
CALL getjd ( jd, year, day, hour )
C
C
time = jd - 51545.0
C
C
mnlong = 280.460 + 0.9856474 * time
mnlong = mod(mnlong,360.0)
IF ( mnlong < 0.0 ) mnlong = mnlong + 360.0
C
C
mnanom = 357.528 + 0.9856003 * time
mnanom = mod(mnanom,360.0)
IF ( mnanom < 0.0 ) mnanom = mnanom + 360.0
mnanom = mnanom * rpd
C
C
eclong = mnlong + 1.915 * SIN mnanom
eclong = eclong + 0.020 * SIN ( 2.0 * mnanom )
eclong = mod(eclong,360.0)
IF ( eclong < 0.0 ) eclong = eclong + 360.0
C
oblqec = 23.439 - 0.0000004 * time
eclong = eclong * rpd
oblqec = oblqec * rpd
C
C
num = COS ( oblqec ) * SIN ( eclong )
den = COS ( eclong )
ra = ATAN ( num / den )
C
IF ( den < 0.0 ) THEN
ra = ra + pi
ELSE IF ( num < 0.0 ) THEN
ra = ra + twopi
END IF
C
C
C
dec = ASIN ( SIN ( oblqec ) * SIN ( eclong ) )
C
C
gmst = 6.697375 + 0.0657098242 * time + hour
gmst = mod(gmst,24.0)
IF ( gmst < 0.0 ) gmst = gmst + 24.0
C
C
lmst = gmst + long / 15.0
lmst = mod(lmst,24.0)
IF ( lmst < 0.0 ) lmst = lmst + 24.0
C
lmst = lmst * 15.0 * rpd
C
C
ha = lmst - ra
C
temp = 0.0 - pi
IF ( ha < temp ) ha = ha + twopi
IF ( ha > pi ) ha = ha - twopi
C
C el = ASIN ( SIN ( dec ) * SIN ( lat * rpd )+
C COS ( dec ) * COS ( lat * rpd ) * COS ( ha ) )
C
temp = COS ( dec ) * COS ( lat * rpd ) * COS ( ha )
temp = temp + SIN ( dec ) * SIN ( lat * rpd )
el = ASIN ( temp )
C
az = 0.0 - SIN ( ha )
az = az * COS ( dec )
az = az / COS ( el )
az = ASIN ( az )
C
C IF ( (SIN(dec)-SIN(el)*SIN(lat*rpd) ) > 0.0 ) THEN
temp = SIN(dec) - SIN (el) * SIN ( lat * rpd )
IF ( temp > 0.0 ) THEN
IF ( SIN ( az ) < 0.0 ) az = az + twopi
ELSE
az = pi - az
END IF
C
C patch for negative azimuths
C probably better to use pi/2 offset for smoothing
C
IF ( az < 0.0 ) az = az + twopi
C
el = el / rpd
az = az / rpd
C
C Refraction as NOAA
C
twil = 0.0 - 0.575
IF ( el > 85.0 ) THEN
refrac = 0.0
ELSE
eltan = TAN ( el * rpd )
IF ( el > 5.0 ) THEN
refrac = 58.1 / eltan
temp = eltan ** 3
temp = 0.07 / temp
refrac = refrac - temp
temp = eltan ** 5
temp = 0.000086 / temp
refrac = refrac + temp
ELSE IF ( el > twil ) THEN
temp = 0.0 - 12.79 + el * 0.711
temp = 103.4 + el * temp
temp = 0.0 - 518.2 + el * temp
refrac = 1735.0 + el * temp
ELSE
refrac = 20.774 / eltan
refrac = 0.0 - refrac
END IF
refrac = refrac / 3600.0
END IF
C
el = el + refrac
C
temp = hour + zone
C
WRITE (*,0) temp, el, az
C
10 CONTINUE
PRINT *, "Finish"
END
C
C
SUBROUTINE getjd ( jd, year, day, hour )
INTEGER*4 delta, leap, jtemp
DOUBLE gtemp
delta = year - 1949
leap = delta / 4
jtemp = delta * 365 + leap + day
gtemp = DBLE ( jtemp )
jd = gtemp + 32916.5 + hour / 24.0
RETURN
END
C
DOUBLE FUNCTION mod ( val, div )
DOUBLE tempval, tempres
tempval = val / div
tempres = DBLE ( IFIX tempval )
tempres = tempres * div
mod = val - tempres
RETURN
END
C