# heliostat program

Discussion in 'Home Power and Microgeneration' started by David Williams, Jul 10, 2007.

1. ### David WilliamsGuest

I recently wrote a QBasic program that calculates the position of the
sun in the sky, as azimuth and elevation, as seen from anywhere on
earth, at any time and date. It also calculates the orientation in which
a mirror should be pointed in order to reflect sunlight in any desired
direction, as is done by a heliostat. With a few hardware-specific
additions to make the computer control the mirror, this program could
be used to run a heliostat. In fact, it is very similar to the program
that I used to run a heliostat, years ago.

For anyone who may be interested, I'll append the program to this
message.

Incidentally, soon after I wrote it, I posted a version of this program
in another newsgroup. There was nothing wrong with that. It worked just
fine. But I have found a few ways to clean up and simplify the coding,
which are incorporated into this version.

Enjoy!

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

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