'===========================================================================
' Subject: Equation of Time                   Date: 06-14-03 (  :  )       
'  Author: David Williams                     Code: Qbasic, QB, PDS        
'  Origin: david.williams@ablelink.org      Packet: DATETIME.ABC
'===========================================================================
' ETIMSDEC.BAS
' David Williams, 2003
' david.williams@ablelink.org

' This version dated 2003 June 14

' Compares results of (home brewed) Equation of Time
' and Solar Declination calculations (performed in
' functions Equt, ETApprox and Declin) with published values.

' Equation of Time is difference between Solar Time and Mean Time.
' Sundials show solar time. Clocks show mean time.
' Solar Declination is latitude of the sun in celestial coordinates.

' Published values are from the book:
' Sundials, Their Theory and Construction
' by: Albert E. Waugh
' Dover Publications, New York, 1973
' ISBN 0-486-22947-5

DECLARE FUNCTION Equt! (D!)
DECLARE FUNCTION Declin! (D!)
DECLARE FUNCTION ETApprox! (D!)

DIM ML(1 TO 13)
FOR X = 1 TO 13
  READ ML(X)
NEXT
DATA 33,29,33,31,33,31,33,33,31,33,31,33,0
' Month lengths adjusted for numbers of pixels on screen

ETData:

DATA 9, 128, 375, -20, 15: ' plotting parameters for E. of T. graph

DATA 3.12,7.38,11.08,13.33,14.19,13.49,12.34,10.18,7.28
DATA 4.08,1.16,-1,-2.51,-3.4,-3.34,-2.25,-.39,1.28
DATA 3.33,5.16,6.15,6.16,5.14,3.16,.12,-3.08,-6.4
DATA -10.05,-13.02,-15.12,-16.2,-16,-14.16,-11.11,-7.02,-2.13
' Above data show published values of Equation of Time on
' 1st, 11th and 21st days of months. Decimal point separates
' minutes, to left, from seconds, to right, e.g. 3.16 means
' 3 minutes and 16 seconds. Seconds are multiplied by 5/3 to
' get fractional minutes after READ instruction.

DeclnData:

DATA 5, 65, 407, -25, 25: ' plotting parameters for Decl'n graph

DATA -23.04,-21.56,-20.05,-17.2,-14.18,-10.52,-7.49,-3.57,0
DATA 4.18,8.07,11.39,14.54,17.43,20.04,21.58,23.02,23.26
DATA 23.09,22.11,20.36,18.1,15.27,12.19,8.3,4.48,.57
DATA -2.57,-6.48,-10.29,-14.14,-17.15,-19.47,-21.43,-22.57,-23.26
' Above data show published values of solar declination on 1st,
' 11th and 21st days of month. See note after Equation of Time data.

SCREEN 12
CL0 = 7: CL1 = 14: CL2 = 11   ' colours used

Menu:
CLS
COLOR CL0
PRINT "1. Draw Equation of Time Graph (Fairly Accurate Function)"
PRINT "2. Draw Equation of Time Graph (Simple Approximate Function)"
PRINT "3. Draw Solar Declination Graph"
PRINT "4. Quit program"
PRINT
PRINT "Which (by number)?";
WHILE INKEY$ <> "": WEND
DO
  F% = VAL(INKEY$)
LOOP UNTIL F% >= 1 AND F% <= 4
PRINT F%
IF F% = 4 GOTO Wayout

' Set up graph
CLS

IF F% < 3 THEN
  PRINT TAB(30); "EQUATION OF TIME"
  LOCATE 3, 15
  PRINT "Graph shows difference in minutes between clock"
  PRINT TAB(15); "and sundial time. Positive difference means"
  PRINT TAB(15); "clock is ahead of sundial, and vice versa."
  LOCATE 26, 15
  COLOR CL1
  PRINT "Graph shows results of ";
  IF F% = 2 THEN PRINT "(approximate) ";
  PRINT "calculation."
  COLOR CL2
  PRINT TAB(15); "Circles show published values of Equation of Time."
  RESTORE ETData
ELSE
  PRINT TAB(25); "Declination of Sun, in Degrees"
  LOCATE 3, 10
  COLOR CL1
  PRINT "Graph shows calculated function. ";
  COLOR CL2
  PRINT "Circles show published values.";
  RESTORE DeclnData
END IF

READ ML%, VT%, VB%, LL%, UL%

COLOR CL0
LOCATE 16, 67
PRINT "-=";

FOR T = LL% TO UL% STEP 5
  LINE (137, 247 - 6.4 * T)-(530, 247 - 6.4 * T)
  LOCATE 16 - T / 2.5, 14
  PRINT RIGHT$(" " + STR$(T), 3);
  IF T = 0 THEN PRINT " =";
NEXT
LOCATE ML%, 20
PRINT "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec";

' Do graphing
T = 0
D = 0
MC = 0
FOR M = 0 TO 384
  M1 = 147 + M
  IF M = T THEN
    ' Draw vertical line
    COLOR CL0
    LINE (M1, VT%)-(M1, VB%)
    MC = MC + 1
    T = T + ML(MC)
    MS = 0
  ELSE
    ' Plot point(s)
    COLOR CL1
    SELECT CASE F%
      CASE 1: G = Equt(D)
      CASE 2: G = ETApprox(D)
      CASE ELSE: G = Declin(D)
    END SELECT
    PSET (M1, 247 - 6.4 * G)
    IF MS = 0 OR MS = 10 OR MS = 20 THEN
      ' Plot circle representing published value
      COLOR CL2
      READ E    ' published values
      E = FIX(E) + (E - FIX(E)) * 5 / 3
      CIRCLE (M1, 247 - 6.4 * E), 2
    END IF
    MS = MS + 1
    D = D + 1
  END IF
  IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1
NEXT

' Tidy up
COLOR CL0
LOCATE 30, 25
PRINT "*** Press any key to continue ***";
WHILE INKEY$ <> "": WEND
WHILE INKEY$ = "": WEND
GOTO Menu

Wayout:
PRINT
PRINT
PRINT "Note: Using the fairly accurate function, the calculated and published"
PRINT "values of Equation of Time show a Root Mean Square difference of about"
PRINT "3.7 *seconds*. The greatest difference is 6.0 seconds, on July 11."
PRINT "Using the approximate function, the R.M.S. difference is about 10.5"
PRINT "seconds. The largest difference is 19.6 seconds, on March 1."
PRINT
PRINT "The approximate calculation is simpler and faster than the more"
PRINT "accurate one, and is certainly adequate for many purposes, such as"
PRINT "comparing the readings of sundials and clocks."
PRINT
PRINT "Comparison of the calculated and published values of Solar Declination"
PRINT "on same dates shows a Root Mean Square difference of 4.7 arc-minutes."
PRINT "The largest difference is 8.8 arc-minutes, on April 11."
PRINT
PRINT "In order to get even better accuracy than this declination calculation"
PRINT "and the more accurate Equation of Time one, much more elaborate"
PRINT "methods are needed."

END

FUNCTION Declin (D)
  ' Calculates declination of sun, in degrees, on Dth day of year.
  ' D = 0 on January 1.
  PI = 3.14159265#
  W = PI / 182.5 ' mean orbital angular speed in radians per day
  ' Solstice 10 days before New Year's Day. Perihelion 2 days after.
  A = W * (D + 10) + .0334 * SIN(W * (D - 2))  ' .0334 is twice eccentricity
  S = -.3979 * COS(A) ' .3979 is sine of axial tilt
  Declin = ATN(S / SQR(1 - S * S)) * 180 / PI ' arc-sine, in degrees
END FUNCTION

FUNCTION Equt (D)
  ' Calculates equation of time, in minutes, on
  ' Dth day of year. (D = 0 on January 1.)
  PI = 3.14159265#
  W = PI / 182.5  ' mean orbital angular speed in radians per day
  T = D + 10 'Solstice 10 days before Jan 1. Perihelion 12 days after solstice
  A = W * T + .0334 * SIN(W * (T - 12)) ' .0334 is twice eccentricity
  B = ATN(TAN(A) / .9174) - W * T ' .9174 is cosine of axial tilt
  B = B - PI * INT(B / PI + .5) ' ATN uncertain. B must be close to zero.
  Equt = B * 720 / PI ' 720 minutes in 12 hours
END FUNCTION

FUNCTION ETApprox (D)
  ' Calculates approximate equation of time, in minutes, on
  ' Dth day of year. (D = 0 on January 1.)
  W = 3.14159265# / 182.5  ' mean orbital angular speed in radians/day
  ' solstice 10 days before new year's day. perihelion 2 days after
  A = W * (D + 10) + .0334 * SIN(W * (D - 2)) ' .0334 is twice eccentricity
  ETApprox = 9.8 * SIN(A + A) + 7.7 * SIN(A - .21) ' approximate formula
END FUNCTION