Kipp&Zonen BSRN Scientific Solar Monitoring System - Manual (Page 181)

169
HC1 = .0001184#
' Constant for the calculation of airm ass
AC1 = -1.253#
' Get the current julian date (actually add 2,400,000 for JD).
Delta = year - basedate
Leap = INT(Delta / 4)
JulianDy = baseday + Delta * ThreeSixtyFive + Leap + jd + GMT / Twentyfour
' First num ber is m id. 0 jan 1949 m inus 2.4e6; Leap = Leap days since 1949.
' Calculate ecliptic coordinates.
Tim e = JulianDy - FiveOneFiveFourFive
' 51545.0 + 2.4e6 = noon 1 jan 2000.
' Force m ean longitude between 0 and 360 degrees.
MnLon = C1 + C2 * Tim e
MnLon = MnLon - INT(MnLon / ThreeSixty) * ThreeSixty
IF MnLon < 0 THEN MnLon = MnLon + ThreeSixty
' Mean anom aly in radians between 0 and 2*Pi
MnAnom = C3 + C4 * Tim e
MnAnom = MnAnom - INT(MnAnom / ThreeSixty) * ThreeSixty
IF MnAnom < 0 THEN MnAnom = MnAnom + ThreeSixty
MnAnom = MnAnom * ToRad
' Com pute ecliptic longitude and obliquity of ecliptic in radians.
EcLon = MnLon + C5 * SIN(MnAnom ) + Point02 * SIN(Two * MnAnom )
EcLon = EcLon - INT(EcLon / ThreeSixty) * ThreeSixty
IF EcLon < 0 THEN EcLon = EcLon + ThreeSixty
OblqEc = C6 - C7 * Tim e
EcLon = EcLon * ToRad
OblqEc = OblqEc * ToRad
' Calculate right ascension and declination.
Num = COS(OblqEc) * SIN(EcLon)
Den = COS(EcLon)
Ra = ATN(Num / Den)
IF Den < 0 THEN
Ra = Ra + pi
ELSEIF Num < 0 THEN
Ra = Ra + TwoPi
END IF
' Declination in radians.
Dec = SIN(OblqEc) * SIN(EcLon)
Dec = ATN(Dec / SQR(One - Dec * Dec))
' Declination in degrees
Decdegrees = Dec * ToDegrees
' Calculate Greenwich m ean sidereal tim e in hours.
GMST = C8 + C9 * Tim e + GMT
' GMT not changed to sidereal since 'tim e' includes the fractional day.
GMST = GMST - INT(GMST / Twentyfour) * Twentyfour
IF GMST < 0 THEN GMST = GMST + Twentyfour
' Calculate local m ean sidereal tim e in radians.
LMST = GMST - Lon / Fifteen
LMST = LMST - INT(LMST / Twentyfour) * Twentyfour
IF LMST < 0 THEN LMST = LMST + Twentyfour
LMST = LMST * Fifteen * ToRad