BBC BASIC
« Sunrise and sunset calculations »

Welcome Guest. Please Login or Register.
Jan 20th, 2018, 4:16pm


Cross-platform BBC BASIC (Win32, Linux x86, Android, Mac OS-X, Raspberry Pi)

« Previous Topic | Next Topic »
Pages: 1  Notify Send Topic Print
 thread  Author  Topic: Sunrise and sunset calculations  (Read 136 times)
Richard Russell
Administrator
ImageImageImageImageImage


member is offline

Avatar




Homepage PM


Posts: 689
xx Sunrise and sunset calculations
« Thread started on: Jun 21st, 2017, 07:02am »

I published this code on the Winter Solstice in 2014, but that was over at the forum that I no longer have access to. So for the Summer Solstice here it is again (it needs DATELIB to be installed).

Code:
      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 

Richard.
User IP Logged

Pages: 1  Notify Send Topic Print
« Previous Topic | Next Topic »

Donate $6.99 for 50,000 Ad-Free Pageviews!


This forum powered for FREE by Conforums ©
Sign up for your own Free Message Board today!
Terms of Service | Privacy Policy | Conforums Support | Parental Controls