World of Solar Thermal.com provides the latest news, views and reviews of the Solar Thermal Marketplace Worldwide
Homepage
Solar Thermal Events Calendar
Forum
World Renewable Energy Association (WREA) Website
Contact Us
Solar Thermal Courses
Renewable Energy Jobs.net

Go Back   World of Solar Thermal.com - The World's #1 Solar Thermal Energy Site! > Main Category > Main Forum

Reply
 
Thread Tools Display Modes
  #1  
Old 08-19-2008
david w david w is offline
Junior Member
 
Join Date: Aug 2008
Posts: 2
Default basic heliostat and sun tracker program

In the late 1980s, I designed, built and programmed a
computer-controlled heliostat. It worked excellently with almost
no attention for many years. The computer was a Commodore VIC 20,
which was old even then, and had only 4.5 kilobytes of memory. Its
particular program would work only on a VIC. However, I have recently
taken the astronomical and trigonometrical parts of the program
and made them into a new program which I'll append below..

I'll append two versions of the program. The first is in QBasic, and
contains quite a lot of explanatory comments. The second version is in
a very generic BASIC, and has been tested on many implementations of
the language. It even has line numbers!

dow

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

' SunAlign.BAS (Version for QBasic and similar dialects)

' Calculates position of sun in sky, as azimuth (compass bearing
' measured clockwise from True North) 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
' 3605 Lakeshore Blvd. West
' Toronto, Ontario. M8W 4Y6
' Canada

' Initially dated 2007 Jul 07
' This version 2008 Jan 13

' All angles in radians except in i/o routines DegIn and DegOut

DECLARE SUB C2P (X, Y, Z, AZ, EL)
DECLARE SUB P2C (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)
DECLARE SUB DegIn (P$, X)
DECLARE SUB DegOut (P$, X)

CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs
CONST DR = 180 / PY ' degree / radian factor
W = 2 * PY / 365 ' earth's mean orbital angular speed in radians/day
WR = PY / 12' earth's speed of rotation relative to sun (radians/hour)
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
SN = 10 * W ' 10 days from December solstice to New Year (Jan 1)
SP = 12 * W ' 12 days from December solstice to perihelion

CLS

Menu:

PRINT "1. Calculate sun's position"
PRINT "2. Calculate mirror orientation"
PRINT "3. Calculate both"
PRINT "4. Quit program"
PRINT
PRINT "Which? (1 - 4)";

DO
S% = VAL(INKEY$)
LOOP UNTIL S% >= 1 AND S% <= 4
PRINT S%

IF S% = 4 THEN END

' Note: For brevity, no error checks on user inputs

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

PRINT

CL = PY / 2 - LT ' co-latitude

D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of year (D = 0 on Jan 1)

A = W * D + SN ' orbit angle since solstice at mean speed
B = A + E2 * SIN(A - SP) ' angle with correction for eccentricity

C = (A - ATN(TAN(B) / CT)) / PY
SL = PY * (C - INT(C + .5))' solar longitude relative to mean position

C = ST * COS(B)
DC = ATN(C / SQR(1 - C * C)) ' solar declination (latitude)
' arcsine of C. ASN not directly available in QBasic

LD = (HR - TZN + MIN / 60) * WR + SL + LG ' longitude difference

CALL P2C(LD, DC, sX, sY, sZ) ' polar axis (perpend'r to azimuth plane)
CALL C2P(sY, sZ, sX, sAZ, sEL) ' horizontal axis
CALL P2C(sAZ - CL, sEL, sY, sZ, sX) ' rotate by co-latitude

IF sZ < 0 THEN
BEEP
PRINT "Sun Below Horizon"
PRINT
GOTO NewCalc
END IF

IF S% <> 2 THEN ' calculate and display sun's position

CALL C2P(sX, sY, sZ, sAZ, sEL) ' vertical axis

DegOut "Sun's azimuth: ", sAZ
DegOut "Sun's elevation: ", sEL

PRINT

END IF

IF S% > 1 THEN ' calculate and display mirror orientation

PRINT "For target direction of light reflected from mirror:"
DegIn "Azimuth of target direction (degrees)", tAZ
DegIn "Elevation of target direction (degrees)", tEL

PRINT

CALL P2C(tAZ, tEL, tX, tY, tZ) ' target vector X,Y,Z
CALL C2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)
' angle bisection by vector addition

PRINT "Mirror aim direction (perpendicular to surface):"
DegOut "Azimuth: ", mAZ
DegOut "Elevation: ", mEL

PRINT

END IF

NewCalc:

PRINT
PRINT "New Calculation"
PRINT

GOTO Menu


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

SELECT CASE SGN(X)
CASE 1: Ang = ATN(Y / X)
CASE -1: Ang = ATN(Y / X) + PY
CASE ELSE: Ang = SGN(Y) * PY / 2
END SELECT

END FUNCTION

SUB C2P (X, Y, Z, AZ, EL)
' Cartesian to Polar. Convert from X,Y,Z to AZ,EL

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

END SUB

SUB DegIn (P$, X)
' Input angle in degrees and convert to radians

PRINT P$;
INPUT N
X = N / DR

END SUB

SUB DegOut (P$, X)
' converts radians to degrees, rounds to nearest 0.1, and prints

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

END SUB

SUB P2C (AZ, EL, X, Y, Z)
' Polar to Cartesian. Convert from AZ,EL to X,Y,Z

Z = SIN(EL)
C = -COS(EL)
X = C * SIN(AZ)
Y = C * COS(AZ)

END SUB

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

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

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

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

Last edited by david w; 08-20-2008 at 05:56 PM.
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Forum Jump



All times are GMT. The time now is 09:34 AM.
Copyright 2009: Pangea Digital Media Limited; All Rights Reserved
Powered by vBulletin Version 3.6.5
Copyright 2000 - 2009, Jelsoft Enterprises Ltd.