## Fidonet Portal

From: DAVID WILLIAMS (1:250/514)
To: All
Date: Sat, 28.07.07 15:55
heliostat code
I don't seem to have posted this here, although I've done so in some
Usenet newsgroups.

I recently wrote some BASIC code that calculates the position of the
sun in the sky, as azimuth and angle of elevation, as seen from
anywhere on earth, at any time on any date. It also calculates the
direction in which a mirror should be aimed if it is to reflect
sunlight in any desired direction. A heliostat, of course, has a mirror
that does that, and is turned so as to continue reflecting in the right
direction as the sun moves across the sky.

I'll append two versions of the program below. The first is in QBasic,
but will also run in some other implementations. The second is in a
very simple BASIC, and should run on almost any version of the
language. However, it is necessarily messier than the first one.

Enjoy. Let me know if you have any questions.

dow

------------------------------------------------

' SunAlign.BAS
' Calculates position of sun in sky, as azimuth (true
' compass bearing) 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
' M8W 4Y6

' Initially dated 2007 Jul 07
' This version 2007 Jul 22

DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff\$ (X)
DECLARE SUB D2P (X, Y, Z, AZ, EL)
DECLARE SUB P2D (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)

CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs
CONST DR = 180 / PY ' degree / radian factor

CLS
PRINT "Use negative numbers for directions opposite to those shown."
INPUT "Observer's latitude (degrees North)"; LT
INPUT "Observer's longitude (degrees West)"; LG
INPUT "Date (M#,D#)"; Mth%, Day%
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
INPUT "Time Zone (+/- hours from GMT/UT)"; TZN

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination

LD = 15 * (HR - TZN) + (MIN + ET) / 4 - LG 'longitude difference
CALL P2D(LD, DC, sX, sY, sZ)
CALL D2P(sY, sZ, sX, qAZ, qEL)
CALL P2D(qAZ + LT - 90, qEL, sY, sZ, sX)
CALL D2P(sX, sY, sZ, sAZ, sEL)

PRINT
IF sEL < 0 THEN PRINT "Sun Below Horizon": END
PRINT "Sun's azimuth: "; ROff\$(sAZ); " degrees"
PRINT "Sun's elevation: "; ROff\$(sEL); " degrees"

PRINT
PRINT "Calculate heliostat mirror alignment? (y/n) ";
DO
K\$ = UCASE\$(INKEY\$)
LOOP UNTIL K\$ = "Y" OR K\$ = "N"
PRINT K\$

IF K\$ = "Y" THEN
PRINT
INPUT "Azimuth of target direction (degrees)"; tAZ
INPUT "Elevation of target direction (degrees)"; tEL
CALL P2D(tAZ, tEL, tX, tY, tZ)
CALL D2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)

PRINT
PRINT "Mirror aim direction (perpendicular to surface):"
PRINT "Azimuth: "; ROff\$(mAZ); " degrees"
PRINT "Elevation: "; ROff\$(mEL); " degrees"
END IF

END

FUNCTION Ang (X, Y)
' calculates angle between positive X axis and vector to (X,Y)

IF X = 0 THEN
A = SGN(Y) * PY / 2
ELSE
A = ATN(Y / X)
IF X < 0 THEN A = A + PY
END IF
Ang = A

END FUNCTION

SUB D2P (X, Y, Z, AZ, EL)
' convert from X,Y,Z to AZ, EL
' Outputs in degrees

EL = Ang(SQR(X * X + Y * Y), Z) * DR
A = Ang(Y, X) * DR
IF A < 180 THEN AZ = A + 180 ELSE AZ = A - 180

END SUB

FUNCTION ET.Dec (D, F%) STATIC
' Calculates equation of time, in minutes, or solar declination,
' in degrees, on day number D of year. (D = 0 on January 1. ' F% selects function: True (non-zero) for Equation of Time,
' False (zero) for Declination.
' STATIC means variables are preserved between calls of function
' This version assumes PY (PI) and DR (180/PI) are initialized

IF W = 0 THEN ' first call, initialize constants

W = 2 * PY / 365 ' earth's mean orbital angular speed in radians/day
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
SP = 12 * W ' 12 days from December solstice to perihelion
D1 = -1 ' holds last D. Saves time if D repeated for both functions

END IF

IF D D1 THEN ' new value of D
A = W * (D + 10) ' Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
D1 = D
END IF

IF F% THEN ' equation of time calculation
C = (A - ATN(TAN(B) / CT) / PY
ET.Dec = 720 * (C - INT(C + .5) ' in 720 minutes, earth rotates PI radians relative to sun

ELSE ' declination calculation
C = ST * COS(B)
ET.Dec = ATN(C / SQR(1 - C * C) * DR
' arcsine of C in degrees. ASN not directly available in QBasic

END IF

END FUNCTION

SUB P2D (AZ, EL, X, Y, Z)
' convert from AZ,EL to X,Y,Z
' Inputs in degrees

E = EL / DR
A = AZ / DR
Z = SIN(E)
C = -COS(E)
X = C * SIN(A)
Y = C * COS(A)

END SUB

FUNCTION ROff\$ (X)
' neatly rounds number to one place of decimals

S\$ = LTRIM\$(STR\$(INT(10 * ABS(X) + .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\$
ROff\$ = LEFT\$(S\$, LEN(S\$) - 1) + "." + RIGHT\$(S\$, 1)

END FUNCTION

----------------------------------------------------

100 REM SunAlign.BAS

110 REM Calculates position of sun in sky, as azimuth (true
120 REM compass bearing) and angle of elevation, as seen from
130 REM 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
180 REM M8W 4Y6

190 REM Initially dated 2007 Jul 07
200 REM This version 2007 Jul 22

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 West)"; 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) / 10: REM roundoff to 1 decimal place
880 RETURN

---------------------------------------------------
--- Platinum Xpress/Win/WINServer v3.0pr5
* Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)