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