REM Based on a program by Roger W. Sinnott
      REM lat = latitude (degrees north)
      REM long = longitude (degrees east)
      REM Z = time zone (hours offset from GMT)
      REM Y = year
      REM M = month (1-12)
      REM D = day (1-31)
      REM tr = sunrise time (hours), -1 indicates no sunrise
      REM ts = sunset time (hours), -1 indicates no sunset
      REM ar = sunrise azimuth (degrees)
      REM as = sunset azimuth (degrees)
      DEF PROCsunriseset(lat, long, Z, Y, M, D, RETURN tr, RETURN ts, RETURN ar, RETURN as)
      LOCAL A0, A1, C, D0, D1, D2, DD, H0, H2, HD, lst, N, S, T, T0, Z0

      T = FN_mjd(D,M,Y) - (FN_mjd(1,1,2000) + 0.5) : REM Days from noon on 1st Jan 2000
      T0 = T / 36525                   : REM Centuries from noon on 1st Jan 2000
      Z0 = -Z / 24                     : REM Time zone offset in days

      lst = 24110.5 + 8640184.812999999 * T0 + 86636.6 * Z0 + 240 * long : REM seconds
      lst = RAD((lst MOD 86400) / 240) : REM Local sidereal time (radians)

      PROCsun(T + Z0, A0, D0)          : REM Sun's position 'today'
      PROCsun(T + Z0 + 1, A1, D1)      : REM Sun's position 'tomorrow'
      IF A1 < A0 THEN A1 += PI * 2

      H0 = lst - A0
      HD = 15 * RAD(1.0027379) - (A1 - A0) / 24
      DD = (D1 - D0) / 24
      S = SIN(RAD(lat))
      C = COS(RAD(lat))

      rise = -1
      set = -1
      FOR N = 0 TO 23
        D2 = D0 + DD : H2 = H0 + HD
        PROCtestevent(C, S, N, H0, H2, D0, D2, tr, ts, ar, as)
        D0 = D2 : H0 = H2
      NEXT
      ENDPROC

      REM Test an hour window for an event
      REM H0,H2 = start and end Hour angle
      REM D0,D2 = start and end Declination
      DEF PROCtestevent(C, S, N, H0, H2, D0, D2, RETURN tr, RETURN ts, RETURN ar, RETURN as)
      LOCAL A, B, D, E, H1, D1, R0, R1, R2, U, V, W, Z
      Z = COS(RAD(90.833)) : REM Zenith distance
      H1 = (H2 + H0) / 2   : REM Hour angle at half-hour
      D1 = (D2 + D0) / 2   : REM declination at half-hour
      R0 = S * SIN(D0) + C * COS(D0) * COS(H0) - Z
      R1 = S * SIN(D1) + C * COS(D1) * COS(H1) - Z
      R2 = S * SIN(D2) + C * COS(D2) * COS(H2) - Z
      IF SGN(R0) = SGN(R2) THEN ENDPROC : REM No sunrise or sunset this hour
      A = 2 * R2 - 4 * R1 + 2 * R0
      B = 4 * R1 - 3 * R0 - R2
      D = B * B - 4 * A * R0 : IF D < 0 THEN ENDPROC : REM No solution to quadratic
      D = SQR(D)
      E = (0 - B + D) / (2 * A)
      IF E > 1 OR E < 0 THEN E = (0 - B - D) / (2 * A)
      U = H0 + E * (H2 - H0)
      V = 0 - COS(D1) * SIN(U)
      W = C * SIN(D1) - S * COS(D1) * COS(U)
      Z = DEG(ATN(V / W))
      IF W < 0 THEN Z += 180
      IF Z < 0 THEN Z += 360
      IF Z > 360 THEN Z -= 360
      IF R0 < 0 THEN tr = N + E : ar = Z ELSE ts = N + E : as = Z
      ENDPROC

      REM Sun's position in the sky (Van Flandern & Pulkkinen, 1979)
      REM D = Delta days from noon on 1st January 2000
      DEF PROCsun(D, RETURN RA, RETURN Dec)
      LOCAL L, G, T, U, V, W
      T = D / 36525 + 1 : REM Number of centuries since 1900
      L = 0.779072 + 0.00273790931 * D
      G = 0.993126 + 0.00273777850 * D
      L = (L - INT(L)) * 2 * PI : REM Mean latitude (radians)
      G = (G - INT(G)) * 2 * PI : REM Mean anomaly (radians)
      U = 1 - 0.03349 * COS(G) - 0.00014 * COS(2 * L) + 0.00008 * COS(L)
      V = (0.39785 - 0.00021 * T) * SIN(L) - 0.01 * SIN(L - G) + 0.00333 * SIN(L + G)
      W = (0.03211 - 0.00008 * T) * SIN(G) - 0.04129 * SIN(2 * L) + \
      \   0.00104 * SIN(2 * L - G) - 0.00035 * SIN(2 * L + G) - 0.0001
      REM Compute Sun's Right Ascension and Declination (radians)
      RA = L + ASN(W / SQR(U - V * V))
      Dec = ASN(V / SQR(U))
      ENDPROC 

INSTALL @lib$+"DATELIB"