## Fidonet Portal

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

' Toronto, Canada

' 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)

IF sEL < 0 THEN PRINT "Sun Below Horizon": END

PRINT "Sun's azimuth: "; ROff$(sAZ); " degrees"

PRINT "Sun's elevation: "; ROff$(sEL); " degrees"

PRINT "Calculate heliostat mirror alignment? (y/n) ";

DO

K$ = UCASE$(INKEY$)

LOOP UNTIL K$ = "Y" OR K$ = "N"

PRINT K$

IF K$ = "Y" THEN

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 "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

170 REM Toronto, Canada

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)