CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
PROGRAM gpssolpos
C
C Reads GPS record and calculates sun position
C
C based on "sunae" subroutine
C (Michalski, J./Harrison, L./mod Wiscombe, W )
C
C Input :
C GPS NMEA $GPRMC string
C
C
C Some factors used for the calculation :
C Time Hours in decimal eg 13.75
C Day,month,year Date eg 1,12,2015
C = 1 Dec 2015 (range 2000 - 2050)
C latitude, longitude S = -ve, E = +ve
C Time zone offset Hours ahead of GMT
C eg Sydney is 10.0.
C
C Output :
C
C Hr Local time
C
C El Elevation ( degrees up from horizon )
C Az Azimuth ( degrees cw from north )
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DOUBLE timedb, lat, lon, x, y
DOUBLE az, el, zone, hourdp, jd
INTEGER*4 timeint, hournt, rmin, min, sec
INTEGER*4 date, year, mnth, day, loclhr
INTEGER*4 latdeg, latmin, latsec
INTEGER*4 londeg, lonmin, lonsec
INTEGER*2 dy, yearday
CHARACTER*20 a, b, c, d, f, h, p, j
CHARACTER*1 e, g
C
C
C
PRINT *, "send GPMRC record now>"
READ a, b, c, d, e, f, g, h, p, j
IF ( a != '$GPMRC' ) THEN
PRINT *, "GPS record not found"
ELSE IF ( c != 'A' ) THEN
PRINT *, "Data invalid"
ELSE
timedb = ATODB b
timeint = IFIX timedb
hournt = timeint / 10000
rmin = timeint - hournt * 10000
min = rmin / 100
sec = rmin - min * 100
C
date = ATOI j
day = date / 10000
mnth = ( date - day * 10000 ) / 100
year = ( date - day * 10000 ) - ( mnth * 100 )
year = year + 2000
C
lat = ATODB d
lon = ATODB f
C
latdeg = IFIX lat / 100
x = lat - DBLE ( latdeg * 100 )
x = x / 60.0
x = DBLE ( latdeg ) + x
IF ( e = 'S' ) THEN
x = 0.0 - x
END IF
C
londeg = IFIX lon / 100
y = lon - DBLE ( londeg * 100 )
y = y / 60.0
y = DBLE ( londeg ) + y
IF ( g = 'W' ) THEN
y = 0.0 - y
END IF
C
lat = x
lon = y
C
PRINT *, "Time (GMT) = ", hournt, ":", min, ":", sec
PRINT *, "Date = ", day, ":", mnth, ":",year
PRINT *, "Lat = ", lat
PRINT *, "Lon = ", lon
END IF
PRINT *, "Time Elev Azim"
C
hourdp = DBLE hournt + DBLE min / 60.0 + DBLE sec / 3600.0
C
CALL sunpos ( hourdp, day, mnth, year, lat, lon, az, el )
C
zone = 10.0
loclhr = hournt + IFIX zone
C
FORMAT (I02,A,I02,A,I02,X,F9.4,X,F9.4,/)
WRITE (*,0) loclhr, ":", min, ":", sec, el, az
C
PRINT *, "Finish"
END
C
C
C
SUBROUTINE sunpos ( hour, dy, mnth, year, lat, long, az, el )
DOUBLE jd, dec, den, eclong, gmst, ha, lmst
DOUBLE mnanom, mnlong, num, oblqec, pi, ra
DOUBLE rpd, refrac, time, twopi, temp, eltan, twil
INTEGER*4 i, day, jtemp
INTEGER*4 delta, leap
INTEGER*2 yr, rem
INTEGER*2 leapday, yearday
C
pi = 3.14159265359
twopi = 2.0 * pi
rpd = 0.01745329252
C
CALL getyearday ( dy, mnth, year, yearday )
CALL getjd ( jd, year, yearday, hour )
C
time = jd - 51545.0
C
mnlong = 280.460 + 0.9856474 * time
mnlong = mod(mnlong,360.0)
IF ( mnlong < 0.0 ) mnlong = mnlong + 360.0
C
mnanom = 357.528 + 0.9856003 * time
mnanom = mod(mnanom,360.0)
IF ( mnanom < 0.0 ) mnanom = mnanom + 360.0
mnanom = mnanom * rpd
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
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
dec = ASIN ( SIN ( oblqec ) * SIN ( eclong ) )
C
gmst = 6.697375 + 0.0657098242 * time + hour
gmst = mod(gmst,24.0)
IF ( gmst < 0.0 ) gmst = gmst + 24.0
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
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 per 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
RETURN
END
C
C
SUBROUTINE getyearday ( dy, mnth, year, yearday )
INTEGER*2 rem, leapday
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 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
leapday = 0
rem = year - ( year / 4 ) * 4
IF ( rem = 0 ) leapday = 1
IF ( mnth > 2 ) yearday = yearday + leapday
C PRINT *, "Day of year is ", yearday
RETURN
END
C
C
SUBROUTINE getjd ( jd, year, yearday, hour )
INTEGER*4 delta, leap, jtemp
DOUBLE gtemp
delta = year - 1949
leap = delta / 4
jtemp = delta * 365 + leap + yearday
gtemp = DBLE ( jtemp )
jd = gtemp + 32916.5 + hour / 24.0
RETURN
END
C
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
\

send GPMRC record now>
$GPRMC,062400.000,A,3130.0000,S,15100.0000,E,0.19,115.89,011215,,,A*70
Time (GMT) = 6:24:0
Date = 1:12:2015
Lat = -31.50000000
Lon = 151.00000000
Time Elev Azim
16:24:00 27.8816 260.3376
Finish