loading

heliostat and sun tracker basic program


In the late 1980s, I designed, built and programmed a
computer-controlled heliostat. Its mirror reflected sunlight into my
living-room, making it much brighter. It worked excellently with almost
no attention for many years, until a neighbour's tree grew and blocked
sunlight from reaching the mirror. The computer was a Commodore VIC 20,
which was old even then, and had only 4.5 kilobytes of memory. The
program I wrote, in Commodore BASIC, fitted into that space and handled
all the control functions. It even included a few "bells and whistles".
For example, at night-time the mirror was automatically parked face
downward to reduce the buildup of dust.

That particular program would work only on a VIC, and I haven't seen
any of those for many years. However, I have recently taken the
astronomical and trigonometrical parts of the program and made them
into a new program which I'll append below. It calculates the position
of the sun in the sky, as azimuth (true compass bearing) and angle of
elevation, as seen from anywhere on the earth at any time on any date.
It also calculates the required orientation of a mirror if it is to
reflect sunlight in any desired direction. With the addition of some
code to enable the computer to control motors, this could become the
software for a computerized sun-tracker or heliostat.

I'll append two versions of the program. The first is in QBasic, and
contains quite a lot of explanatory comments. The second version is in
a very generic BASIC, and has been tested on many implementations of
the language. It even has line numbers! Personally, I prefer the QBasic
version. The coding is more elegant. However, the generic version is
likely to be useful to more people.

It's public-domain. Use it for any purpose, even commercially.

Enjoy!

dow



' SunAlign.BAS (Version for QBasic and similar dialects)

' Calculates position of sun in sky, as azimuth (compass bearing
' measured clockwise from True North) and angle of elevation, as
' seen from any place on earth, on any date and any time.
' Also calculates alignment of a heliostat mirror.

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

' Initially dated 2007 Jul 07
' This version 2008 Jan 13

' All angles in radians except in i/o routines DegIn and DegOut

DECLARE SUB C2P (X, Y, Z, AZ, EL)
DECLARE SUB P2C (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)
DECLARE SUB DegIn (P$, X)
DECLARE SUB DegOut (P$, X)

CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs
CONST DR = 180 / PY ' degree / radian factor
W = 2 * PY / 365 ' earth's mean orbital angular speed in radians/day
WR = PY / 12' earth's speed of rotation relative to sun (radians/hour)
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
SN = 10 * W ' 10 days from December solstice to New Year (Jan 1)
SP = 12 * W ' 12 days from December solstice to perihelion

CLS

Menu:

PRINT "1. Calculate sun's position"
PRINT "2. Calculate mirror orientation"
PRINT "3. Calculate both"
PRINT "4. Quit program"
PRINT
PRINT "Which? (1 - 4)";

DO
S% = VAL(INKEY$)
LOOP UNTIL S% >= 1 AND S% <= 4
PRINT S%

IF S% = 4 THEN END

' Note: For brevity, no error checks on user inputs

PRINT
PRINT "Use negative numbers for directions opposite to those shown."
PRINT
DegIn "Observer's latitude (degrees North)", LT
DegIn "Observer's longitude (degrees East)", LG
INPUT "Time Zone (+/- hours from GMT/UT)"; TZN
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
INPUT "Date (M#,D#)"; Mth%, Day%

PRINT

CL = PY / 2 - LT ' co-latitude

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of year (D = 0 on Jan 1)

A = W * D + SN ' orbit angle since solstice at mean speed
B = A + E2 * SIN(A - SP) ' angle with correction for eccentricity

C = (A - ATN(TAN(B) / CT)) / PY
SL = PY * (C - INT(C + .5))' solar longitude relative to mean position

C = ST * COS(B)
DC = ATN(C / SQR(1 - C * C)) ' solar declination (latitude)
' arcsine of C. ASN not directly available in QBasic

LD = (HR - TZN + MIN / 60) * WR + SL + LG ' longitude difference

CALL P2C(LD, DC, sX, sY, sZ) ' polar axis (perpend'r to azimuth plane)
CALL C2P(sY, sZ, sX, sAZ, sEL) ' horizontal axis
CALL P2C(sAZ - CL, sEL, sY, sZ, sX) ' rotate by co-latitude

IF sZ < 0 THEN
BEEP
PRINT "Sun Below Horizon"
PRINT
GOTO NewCalc
END IF

IF S% <> 2 THEN ' calculate and display sun's position

CALL C2P(sX, sY, sZ, sAZ, sEL) ' vertical axis

DegOut "Sun's azimuth: ", sAZ
DegOut "Sun's elevation: ", sEL

PRINT

END IF

IF S% > 1 THEN ' calculate and display mirror orientation

PRINT "For target direction of light reflected from mirror:"
DegIn "Azimuth of target direction (degrees)", tAZ
DegIn "Elevation of target direction (degrees)", tEL

PRINT

CALL P2C(tAZ, tEL, tX, tY, tZ) ' target vector X,Y,Z
CALL C2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)
' angle bisection by vector addition

PRINT "Mirror aim direction (perpendicular to surface):"
DegOut "Azimuth: ", mAZ
DegOut "Elevation: ", mEL

PRINT

END IF

NewCalc:

PRINT
PRINT "New Calculation"
PRINT

GOTO Menu

FUNCTION Ang (X, Y)
' calculates angle from positive X axis to vector to (X,Y)

SELECT CASE SGN(X)
CASE 1: Ang = ATN(Y / X)
CASE -1: Ang = ATN(Y / X) + PY
CASE ELSE: Ang = SGN(Y) * PY / 2
END SELECT

END FUNCTION

SUB C2P (X, Y, Z, AZ, EL)
' Cartesian to Polar. Convert from X,Y,Z to AZ,EL

EL = Ang(SQR(X * X + Y * Y), Z)
A = Ang(Y, X)
IF A < PY THEN AZ = A + PY ELSE AZ = A - PY

END SUB

SUB DegIn (P$, X)
' Input angle in degrees and convert to radians

PRINT P$;
INPUT N
X = N / DR

END SUB

SUB DegOut (P$, X)
' converts radians to degrees, rounds to nearest 0.1, and prints

S$ = LTRIM$(STR$(INT(10 * ABS(X * DR) + .5)))
IF S$ = "3600" THEN S$ = "0"
IF LEN(S$) = 1 THEN S$ = "0" + S$
IF X < 0 THEN IF VAL(S$) THEN S$ = "-" + S$
PRINT P$; LEFT$(S$, LEN(S$) - 1); "."; RIGHT$(S$, 1); " degrees"

END SUB

SUB P2C (AZ, EL, X, Y, Z)
' Polar to Cartesian. Convert from AZ,EL to X,Y,Z

Z = SIN(EL)
C = -COS(EL)
X = C * SIN(AZ)
Y = C * COS(AZ)

END SUB





100 REM SunAlign.BAS (Generic BASIC version)

110 REM Calculates position of sun in sky, as azimuth (compass bearing
120 REM measured clockwise from True North) and angle of elevation, as
130 REM seen from any place on earth, on any date and any time.
140 REM Also calculates alignment of a heliostat mirror.

150 REM David Williams
160 REM P.O. Box 48512
170 REM 3605 Lakeshore Blvd. West
180 REM Toronto, Ontario. M8W 4Y6
190 REM Canada

200 REM Original date 2007 Jul 07. This version 2007 Oct 07

210 REM Note: For brevity, no error checks on user inputs

220 CLS
230 PRINT "Use negative numbers for opposite directions."
240 INPUT "Observer's latitude (degrees North)"; LT
250 INPUT "Observer's longitude (degrees East)"; LG
260 INPUT "Date (M#,D#)"; Mth, Day
270 INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
280 INPUT "Time Zone (+/- hours from GMT/UT)"; TZN

290 PY = 4 * ATN(1): REM "PI" not assignable in some BASICs
300 DR = 180 / PY: REM degree/radian factor
310 W = 2 * PY / 365: REM earth's mean orbital speed in radians/day
320 C = -23.45 / DR: REM reverse angle of axial tilt in radians
330 ST = SIN(C): REM sine of reverse tilt
340 CT = COS(C): REM cosine of reverse tilt
350 E2 = 2 * .0167: REM twice earth's orbital eccentricity
360 SP = 12 * W: REM 12 days from December solstice to perihelion

370 D = INT(30.6 * ((Mth + 9) MOD 12) + 58.5 + Day) MOD 365
380 A = W * (D + 10): REM Solstice 10 days before Jan 1
390 B = A + E2 * SIN(A - SP)

400 C = (A - ATN(TAN(B) / CT)) / PY
410 ET = 720 * (C - INT(C + .5)): REM equation of time
420 REM in 720 minutes, earth rotates PI radians relative to sun

430 C = ST * COS(B)
440 EL = ATN(C / SQR(1 - C * C)) * DR: REM solar declination

450 AZ = 15 * (HR - TZN) + (MIN + ET) / 4 + LG: REM longitude diff
460 GOSUB 800
470 R = SQR(Y * Y + Z * Z)
480 AX = Y: AY = Z: GOSUB 710
490 A = AA + (90 - LT) / DR
500 Y = R * COS(A)
510 Z = R * SIN(A)
520 GOSUB 740

530 PRINT : REM AZ & EL are now sun's azimuth & elevation in degrees
540 IF EL < 0 THEN PRINT "Sun Below Horizon": END
550 R = AZ: GOSUB 870: PRINT "Sun's azimuth: "; R; " degrees"
560 R = EL: GOSUB 870: PRINT "Sun's elevation: "; R; " degrees"

570 PRINT
580 INPUT "Calculate heliostat mirror alignment (y/n)"; K$
590 IF K$ = "N" OR K$ = "n" THEN END

600 SX = X: SY = Y: SZ = Z

610 PRINT
620 INPUT "Azimuth of target direction (degrees)"; AZ
630 INPUT "Elevation of target direction (degrees)"; EL
640 GOSUB 800
650 X = X + SX: Y = Y + SY: Z = Z + SZ: GOSUB 740

660 PRINT : REM AZ & EL are now aim azimuth & elevation in degrees
670 PRINT "Mirror aim direction (perpendicular to surface):"
680 R = AZ: GOSUB 870: PRINT "Azimuth: "; R; " degrees"
690 R = EL: GOSUB 870: PRINT "Elevation: "; R; " degrees"

700 END

710 IF AX = 0 THEN AA = SGN(AY) * PY / 2: RETURN
720 AA = ATN(AY / AX): IF AX < 0 THEN AA = AA + PY
730 RETURN

740 AX = SQR(X * X + Y * Y): AY = Z: GOSUB 710
750 EL = AA * DR
760 AX = Y: AY = X: GOSUB 710
770 AZ = AA * DR
780 IF AZ < 180 THEN AZ = AZ + 180 ELSE AZ = AZ - 180
790 RETURN

800 E = EL / DR
810 A = AZ / DR
820 Z = SIN(E)
830 C = 0 - COS(E): REM Won't work without "0" in Liberty Basic
840 X = C * SIN(A)
850 Y = C * COS(A)
860 RETURN

870 R = INT(10 * R + .5): IF R = 3600 THEN R = 0
880 R = R / 10
890 RETURN