'===========================================================================
' Subject: SEASONS CALCULATOR                Date: 10-29-99 (22:14)       
' Author:  Wayne Henderson                   Code: QB, QBasic, PDS        
' Origin:  whenders@becon.org              Packet: DATETIME.ABC
'===========================================================================
'
'   SEASONS CALCULATOR (approximate) times within ñ minute for years
'                                      from +1000 to +3000 or more
'
' from algorithms presented in 'Astronomical Algorithms' by
' Jean Meeus, 1991  ISBN 0-943396-35-2 (get it at the library to see how
'                                       this complicated algorithm works.)
' More accurate methods exist and are mentioned in the book's bibliography
'
DEFINT A-Z
DECLARE FUNCTION dcos# (deg#) 'convert degrees to radians, get cos
DECLARE FUNCTION jd2date$ (month, day AS DOUBLE, y&) 'Julian Day to Calendar date
DIM yr AS DOUBLE, day AS DOUBLE, jde AS DOUBLE, Lambda AS DOUBLE
DIM March AS DOUBLE, June AS DOUBLE, September AS DOUBLE, December AS DOUBLE
DIM c AS INTEGER, b(72) AS DOUBLE
CONST rad# = 3.141592653589793# / 180
'
start: PRINT : INPUT "Year ", q$: yr = VAL(q$): IF yr = 0 THEN SYSTEM
yr = (yr - 2000) / 1000: PRINT
'
' step 1  calculate:
'
'         JDE0 - 2451545.0
'    T =  ----------------, # of Julian centuries since Jan. 1, 2000 noon UT
'               36525
'
'    W = 35999.373T - 2.47 (degrees)
'
'    Lambda = 1 + 0.0334 cos W + 0.0007 cos 2W
'
'  jde0 values
March = 2451623.80984# + 365242.37404# * yr + .05169 * yr * yr
March = March - .00411 * yr * yr * yr - .00057 * yr * yr * yr * yr
June = 2451716.56767# + 365241.62603# * yr + .00325 * yr * yr
June = June + .00888 * yr * yr * yr - .0003 * yr * yr * yr * yr
September = 2451810.21715# + 365242.01767# * yr - .11575 * yr * yr
September = September + .00337 * yr * yr * yr + .00078 * yr * yr * yr * yr
December = 2451900.05952# + 365242.74049# * yr - .06223 * yr * yr
December = December - .00823 * yr * yr * yr + .00032 * yr * yr * yr * yr
'
' step 2  get the 24 terms of the series A cos (B + CT)
'         and calculate their sum in subroutine kalk.
'
RESTORE: FOR c = 1 TO 72: READ a$: b(c) = VAL(a$): NEXT
'  terms of the series:  A cos (B + CT),  S = ä { A cos (B + CT) }
'                 eg: 485 * cos (324.96 + 1934.136 * T) + 203 * ...
DATA 485,324.96,1934.136,203,337.23,32964.467,199,342.08,20.186
DATA 182,27.85,445267.112,156,73.14,45036.886,136,171.52,22518.443
DATA 77,222.54,65928.934,74,296.72,3034.906,70,243.58,9037.513
DATA 58,119.81,33718.147,52,297.17,150.678,50,21.02,2281.226
DATA 45,247.54,29929.562,44,325.15,31555.956,29,60.93,4443.417
DATA 18,155.12,67555.328,17,288.79,4562.452,16,198.04,62894.029
DATA 14,199.76,31436.921,12,95.39,14577.848,12,287.11,31931.756
DATA 12,320.81,34777.259,9,227.73,1222.114,8,15.45,16859.074
'
'........ Change TZ and TZ1$, TZ2$ to the correct values for your time zone:
'........ 8 for PST 7 for PDT, 7 for MST 6 for MDT, 6 for CST 5 for CDT, etc.
'
tz1$ = " EST": tz2$ = " EDT"
jde = March: tz = 5: GOSUB kalk: PRINT "Spring: ";
GOSUB jddate: PRINT jd2date$(month, day, y&); tz1$
jde = June: tz = 4: GOSUB kalk: PRINT "Summer: ";
GOSUB jddate: PRINT jd2date$(month, day, y&); tz2$
jde = September: tz = 4: GOSUB kalk: PRINT "Autumn: ";
GOSUB jddate: PRINT jd2date$(month, day, y&); tz2$
jde = December: tz = 5: GOSUB kalk: PRINT "Winter: ";
GOSUB jddate: PRINT jd2date$(month, day, y&); tz1$
GOTO start
'
' step 3
'                           0.00001S
'             JDE = jde0 + ----------
'                            Lambda
'
kalk:
   t# = (jde - 2451545) / 36525
   w# = 35999.373# * t# - 2.47#
   Lambda = 1 + .0334 * dcos(w#) + .0007 * dcos(2 * w#)
   sum# = 0
   FOR c = 1 TO 72 STEP 3
      sum# = sum# + b(c) * dcos(b(c + 1) + b(c + 2) * t#)
   NEXT
   jde = jde + (.00001 * sum#) / Lambda - tz / 24 'UT to local - tz/24
RETURN
'
jddate: 'algorithm from Peter Baum (pbaum@capecod.net)
   r$ = "312831303130313130313031": s& = INT(jde - 1721118.5#)
   H# = jde - 1721118.5# - s&: m# = 100 * s& - 25
   a = m# \ 3652425: b = a - a \ 4: y& = (100 * b + m#) \ 36525
   c = b + s& - 365 * y& - y& \ 4: month = (5 * c + 456) \ 153
   day = c - (153 * month - 457) \ 5 + H#
   IF month > 12 THEN y& = y& + 1: month = month - 12
   IF y& / 400 = y& \ 400 THEN MID$(r$, 4, 1) = "9"
   IF (y& / 100 <> y& \ 100) AND (y& / 4 = y& \ 4) THEN
      MID$(r$, 4, 1) = "9"
   END IF
RETURN
'
' by Wayne Henderson    Oct. 29, 1999 (Fri)   10:00 pm EDT

DEFSNG A-Z
FUNCTION dcos# (deg#)
   dcos# = COS(deg# * rad#)
END FUNCTION

DEFINT A-Z
FUNCTION jd2date$ (month, day AS DOUBLE, y&)
   r$ = "312831303130313130313031"
   hr! = (day - INT(day)) * 24 'fraction of a day
   min! = (hr! - INT(hr!)) * 60 'fraction of an hour
   sec = CINT((min! - INT(min!)) * 60) 'fraction of a minute
   min! = INT(min!): day = INT(day): hr! = INT(hr!)
   IF sec = 60 THEN sec = 0: min! = min! + 1
   IF min! = 60 THEN min! = 0: hr! = hr! + 1
   IF hr! = 24 THEN hr! = 0: day = day + 1
   IF day > VAL(MID$(r$, month * 2 - 1, 2)) THEN day = 1: month = month + 1
   IF month > 12 THEN month = 1: year = year + 1
   IF month < 10 THEN u$ = "0" ELSE u$ = ""
   IF day < 10 THEN w$ = "0" ELSE w$ = ""
   IF hr! < 10 THEN v$ = "0" ELSE v$ = ""
   IF min! < 10 THEN x$ = "0" ELSE x$ = ""
   IF sec < 10 THEN p$ = "0" ELSE p$ = ""
   u$ = u$ + LTRIM$(STR$(month)) + "-": w$ = w$ + LTRIM$(STR$(day)) + "-"
   v$ = v$ + LTRIM$(STR$(hr!)) + ":": x$ = x$ + LTRIM$(STR$(min!)) + ":"
   p$ = p$ + LTRIM$(STR$(sec))
   jd2date$ = u$ + w$ + LTRIM$(STR$(y&)) + " " + v$ + x$ + p$
END FUNCTION