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

Posted by david williams 10 years ago