#1




basic heliostat and sun tracker program
In the late 1980s, I designed, built and programmed a
computercontrolled heliostat. It worked excellently with almost no attention for many years. The computer was a Commodore VIC 20, which was old even then, and had only 4.5 kilobytes of memory. Its particular program would work only on a VIC. 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.. 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! 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) (24hr format)"; HR, MIN INPUT "Date (M#,D#)"; Mth%, Day% CL = PY / 2  LT ' colatitude 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 colatitude 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) (24hr 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  Last edited by david w; 08202008 at 05:56 PM. 