## 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 "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 "Use negative numbers for directions opposite to those shown."

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%

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"

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

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

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

END IF

NewCalc:

PRINT "New Calculation"

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