As I said, it's in QBasic. It contains some graphics, so it can't easily be completely translated to another dialect, but the FUNCTION coding should be easy to translate.

It's public-domain. Free for any kind of use, including commercial use.

I hope someone finds it useful.

                             dow

---------------------------------------------------------

' ETIMSDEC.BAS
' David Williams, 2003

' P.O. Box 48512
' 3605 Lakeshore Blvd. West
' Toronto, Ontario. M8W 4Y6
' Canada

' This version dated 2007 Jun 13

DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff$ (X)
DECLARE SUB Waitkey (X%)

CL0 = 7: CL1 = 12: CL2 = 15   ' colours used

SCREEN 12

COLOR CL2
CLS

PRINT "ETIMSDEC.BAS"
PRINT
PRINT "Shows Equation of Time and Solar Declination calculations"
PRINT "(performed in function ET.Dec) and compares their results"
PRINT "graphically with published values, showing close agreement."
PRINT
PRINT "These calculations can be used in the programs of computer-"
PRINT "controlled solar-energy equipment, such as sun trackers and"
PRINT "heliostats."
PRINT
PRINT "Equation of Time is difference between Solar Time and Mean"
PRINT "Time. Sundials show solar time. Clocks show mean time. The"
PRINT "Equation of Time is related to solar longitude, relative to"
PRINT "its mean value (the value it would have if the sun's apparent"
PRINT "motion were at uniform speed). One degree westward of solar"
PRINT "longitude corresponds to four minutes in positive direction in"
PRINT "the Equation of Time."
PRINT
PRINT "Solar Declination is latitude of sun in celestial coordinates."
PRINT

' initialize array used in graphing routine
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

DO

  COLOR CL0

  PRINT "1. Find Equation of Time and Solar Declination on given date"
  PRINT "2. Draw Equation of Time Graph"
  PRINT "3. Draw Solar Declination Graph"
  PRINT "4. Draw Solar Analemma"
  PRINT "5. Quit program"
  PRINT
  PRINT "Which (1 - 5)? ";
  WHILE INKEY$ <> "": WEND
  DO
    K$ = INKEY$
  LOOP UNTIL K$ >= "1" AND K$ <= "5"
  PRINT K$

  SELECT CASE VAL(K$)
    CASE 1: GOSUB Calc
    CASE 2, 3: GOSUB Graph
    CASE 4: GOSUB Analemma
    CASE 5: EXIT DO
  END SELECT

LOOP

GOTO WayOut

Calc:

PRINT
COLOR CL1

DO
  INPUT "Date (M#,D#)"; Mth%, Day%
LOOP UNTIL Mth% > 0 AND Mth% < 13 AND Day% > 0 AND Day% < 32

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of 365-day year. Jan 1 = 0. Dec 31 = 364

PRINT
PRINT "Day number"; D + 1; "of year" ' conventional range, 1 - 365

ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination

PRINT
PRINT "Equation of Time: "; ROff$(ET); " minutes"
PRINT "Solar Declination: "; ROff$(DC); " degrees"
PRINT

RETURN ' to menu

Graph:

F% = (K$ = "2") ' true if Equation of Time selected

' Set up graph
CLS

IF F% 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); "sundial is ahead of clock, and vice versa."
  LOCATE 24, 15
  COLOR CL1
  PRINT "Graph shows results of calculation."
  COLOR CL2
  PRINT TAB(15); "Crosses 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 "Crosses 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
    YY = 247 - 6.4 * ET.Dec(D, F%) ' function used
    IF POINT(M1, YY) <> CL2 THEN PSET (M1, YY) ' don't overwrite cross
    IF MS = 0 OR MS = 10 OR MS = 20 THEN
      ' Plot cross representing published value
      COLOR CL2
      READ E ' published value
      E = FIX(E) + (E - FIX(E)) * 5 / 3 ' see note after E of T data
      YY = 247 - 6.4 * E
      LINE (M1 - 2, YY)-(M1 + 2, YY) ' lines make
      LINE (M1, YY - 2)-(M1, YY + 2) ' cross
    END IF
    MS = MS + 1
    D = D + 1
  END IF
  IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1
NEXT

Waitkey 25

RETURN ' to menu

ETData:

DATA 7, 97, 343, -15, 20: ' 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 part is 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 (in degrees and minutes) of solar
' declination on 1st, 11th and 21st days of month. See note after
' Equation of Time data.

' 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

Analemma:

CLS

RESTORE Anadata

PRINT "SOLAR ANALEMMA"
PRINT
PRINT
PRINT "Shows apparent path of sun,"
PRINT "through the year, relative"
PRINT "to its mean position"
PRINT "(marked by cross)."

LOCATE 3, 60
PRINT "N"
LOCATE 7, 60
PRINT "S"
LOCATE 5, 56
PRINT "E       W"

LINE (475, 52)-(475, 90)  ' compass lines
LINE (456, 71)-(494, 71)

LOCATE 9, 1
PRINT TAB(53); "Looking at sky,"
PRINT TAB(53); "if North is upward,"
PRINT TAB(53); "East is to left."

LOCATE 24, 1
PRINT "For dimensions, see graphs"
PRINT "of Solar Declination and"
PRINT "Equation of Time."
PRINT "+4 minutes in Equation of Time"
PRINT "equal 1 degree Westward."

LOCATE 15, 1
PRINT TAB(50); "Celestial"
PRINT TAB(50); "Equator"

X0 = 320  ' origin coordinates
Y0 = 240

LINE (X0, Y0 + 5)-(X0, Y0 - 5) ' cross at (0,0)
LINE (X0 + 5, Y0)-(X0 - 5, Y0)

FOR S = -1 TO 1 STEP 2
  LINE (X0, Y0 + S * 210)-(X0, Y0 + S * 225) ' lines marking N-S
  LINE (X0 + S * 30, Y0)-(X0 + S * 55, Y0)   ' and E-W
NEXT

LOCATE 10, 1
PRINT "Month colours:";

M = 0 ' month number
DM = 0 ' day of month
LM = 0 ' length of month

' draw analemma
FOR DY = 0 TO 364 ' day of year
  IF DM = LM THEN ' start new month
    M = M + 1
    C = M XOR 14 * (M AND 1) ' twiddle bits to get good contrasts
    IF C = 8 THEN C = 14  ' yellow better than dull grey
    COLOR C
    READ LM, MN$
    PRINT TAB(17); MN$ ' month name in correct colour
    DM = 0 ' reset day of month
  END IF
  DM = DM + 1
  FOR D = DY TO DY + .9 STEP .25 ' 4 points per day to avoid gaps
    E = ET.Dec(D, 1) ' equation of time
    L = ET.Dec(D, 0) ' solar declination (latitude)
    PSET (X0 + 2.5 * E, Y0 - 10 * L) ' 4:1 ratio.  4 mins = 1 degree
  NEXT
NEXT

Waitkey 48

RETURN ' to menu

Anadata:

DATA 31,Jan,28,Feb,31,Mar,30,Apr,31,May,30,Jun
DATA 31,Jul,31,Aug,30,Sep,31,Oct,30,Nov,31,Dec
' month lengths and names

WayOut:

COLOR CL2
CLS
PRINT "These calculations of Equation of Time and Solar Declination"
PRINT "are simplified and approximate. However, they are quite good."
PRINT "All differences between calculated and published values of"
PRINT "the solar declination are small compared with the angular"
PRINT "size of the sun in the sky. Similarly, differences between"
PRINT "calculated and published values of the equation of time are"
PRINT "small compared with the time the sun takes to traverse its"
PRINT "own diameter as it moves across the sky. The non-zero size"
PRINT "of the sun, rather than any inaccuracies of calculation,"
PRINT "is the limiting factor on how accurately the sun can be"
PRINT "tracked in solar energy applications, using this routine."
PRINT
PRINT "Leap years cause minor fluctuations in the Equation of Time"
PRINT "and Solar Declination on given dates of the year. However,"
PRINT "the errors that are introduced by ignoring these fluctuations"
PRINT "are always smaller than the apparent radius of the sun. This"
PRINT "program therefore ignores leap years. February 29, when it"
PRINT "occurs, is considered to be the same day as March 1."
PRINT
PRINT "The Equation of Time is treated here as the correction to be"
PRINT "subtracted from a sundial reading to get local mean (clock)"
PRINT "time. It is therefore positive when sundial is ahead of clock."
PRINT "This is the usual sign convention, but the opposite usage is"
PRINT "sometimes found. Take care when comparing values of the"
PRINT "Equation of Time from different sources."

Waitkey 1

PRINT "Note that atmospheric refraction of light affects the apparent"
PRINT "position of the sun in the sky when it is close to the"
PRINT "horizon. Sunlight is then too weak to be used for most"
PRINT "solar-energy applications, so refraction is not usually"
PRINT "included in calculations related to solar energy. This program"
PRINT "does not take account of refraction."
PRINT
PRINT "Comments in the function ET.Dec show how quantities in it are"
PRINT "related to astronomical parameters. Over long periods of time"
PRINT "(many decades or more), some of these parameters will change"
PRINT "significantly. The function should be modified accordingly if"
PRINT "it is to be used over very long periods of time."

Waitkey 1

END

FUNCTION ET.Dec (D, F%) STATIC
  ' Calculates equation of time, in minutes, or solar declination,
  ' in degrees, on day number D of year. (D = 0 on January 1.)
  ' F% selects function: True (non-zero) for Equation of Time,
  ' False (zero) for Declination.
  ' STATIC means variables are preserved between calls of function

IF PI = 0 THEN ' first call, initialize constants

 PI = 4 * ATN(1)
 W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
 DR = 180 / PI ' degree/radian factor
 C = -23.45 / DR ' reverse angle of earth's axial tilt in radians
 ST = SIN(C) ' sine of reverse tilt
 CT = COS(C) ' cosine of reverse tilt
 E2 = 2 * .0167 ' twice earth's orbital eccentricity
 SP = 12 * W ' 12 days from December solstice to perihelion
 D1 = -1 ' holds last D. Saves time if D repeated for both functions

END IF

IF D <> D1 THEN ' new value of D
  A = W * (D + 10) ' Solstice 10 days before Jan 1
  B = A + E2 * SIN(A - SP)
  D1 = D
END IF

IF F% THEN  ' equation of time calculation
  C = (A - ATN(TAN(B) / CT)) / PI
  ET.Dec = 720 * (C - INT(C + .5))
  ' in 720 minutes, earth rotates PI radians relative to sun

ELSE ' declination calculation
  C = ST * COS(B)
  ET.Dec = ATN(C / SQR(1 - C * C)) * DR
  ' arcsine of C in degrees. ASN not directly available in QBasic

END IF

END FUNCTION

FUNCTION ROff$ (X)
 ' neatly rounds number to one place of decimals

S$ = LTRIM$(STR$(INT(10 * ABS(X) + .5)))
IF LEN(S$) = 1 THEN S$ = "0" + S$
IF VAL(S$) THEN S$ = MID$("-x+", SGN(X) + 2, 1) + S$
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1)

END FUNCTION

SUB Waitkey (X%)
 ' prints prompt at position X% on bottom row, then waits

COLOR 15
LOCATE 30, X%
PRINT "* Press any key to continue *";
WHILE INKEY$ <> "": WEND
WHILE INKEY$ = "": WEND
CLS

END SUB