heliostat program

[David Williams]

I wrote this QBasic program a few days ago. It calculates the position
of the sun in the sky as azimuth (true compass bearing) and angle of
elevation, as seen from anywhere on earth, at any time on any date. It
also calculates the orientation of a mirror if it is to reflect
sunlight in any desired direction, as is done in a heliostat. With some
additions to make the computer control motors that turn the mirror,
this program could be the basis of software for a computer-controlled
heliostat.

I hope someone finds it useful.

dow

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

' SunAlign
' 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 09

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 PI = 3.1415926536#
CONST DR = 180 / PI ' 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)
RR = SQR(sY * sY + sZ * sZ)
CL = (90 - LT) / DR ' colatitude in radians
AA = Ang(sY, sZ) + CL
sY = RR * COS(AA)
sZ = RR * SIN(AA)
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) * PI / 2
ELSE
A = ATN(Y / X)
IF X < 0 THEN A = A + PI
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 PI and DR (180/PI) are already initialized

IF W = 0 THEN ' first call, initialize constants

W = 2 * PI / 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)) / PI
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

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

[Morris Dovey]

David Williams wrote:

| I wrote this QBasic program a few days ago. It calculates the
| position of the sun in the sky as azimuth (true compass bearing)
| and angle of elevation, as seen from anywhere on earth, at any time
| on any date. It also calculates the orientation of a mirror if it
| is to reflect sunlight in any desired direction, as is done in a
| heliostat. With some additions to make the computer control motors
| that turn the mirror, this program could be the basis of software
| for a computer-controlled heliostat.

<code snipped>

What have you changed/corrected from the previous version?

--
Morris Dovey
DeSoto Solar
DeSoto, Iowa USA
http://www.iedu.com/DeSoto/

[David Williams]

-> What have you changed/corrected from the previous version?

There are no corrections. As I said, the initial version worked fine.
All I did was to streamline things in some places. For example, the
initial version had a FUNCTION called Ang2, which was called only once.
So I just wrote the code directly into the routine that called it, and
deleted the function. There were a few other simplifications, too.

If you compare the two versions, you'll find the changes. But if you've
used the first version for something and are wondering whether you must
switch, don't worry about it unless you prefer neat code.

dow