heliostat and sun tracker basic program Answered
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.
' 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
' 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
PRINT "1. Calculate sun's position"
PRINT "2. Calculate mirror orientation"
PRINT "3. Calculate both"
PRINT "4. Quit program"
PRINT "Which? (1 - 4)";
S% = VAL(INKEY$)
LOOP UNTIL S% >= 1 AND S% <= 4
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
PRINT "Sun Below Horizon"
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
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
PRINT "New Calculation"
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
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
SUB DegIn (P$, X)
' Input angle in degrees and convert to radians
X = N / DR
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"
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)
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
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"
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
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"
710 IF AX = 0 THEN AA = SGN(AY) * PY / 2: RETURN
720 AA = ATN(AY / AX): IF AX < 0 THEN AA = AA + PY
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
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)
870 R = INT(10 * R + .5): IF R = 3600 THEN R = 0
880 R = R / 10