Here's a program that should be useful in the context of solar energy. It does the astronomical and trigonometrical calculations that have to be done in a computer-controlled sun-tracker or heliostat. With the addition of code to allow the computer to control motors that turn the tracker or mirror, the program could become the software for such a machine.
I'll post two versions. The first is in QBasic but will also run under Borland Turbo BASIC and some other implementations. The second is in a very generic BASIC, with line numbers. It will run (and has been tested) under many varieties of BASIC. The QBasic version is the "real" program, which should be easier to read and understand than the other.
If anyone wants to translate it to another language and post the result, please go ahead!
This is public-domain, free for any kind of use, including commercial use.
Enjoy!
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
--------------------------------------------------------------
--------------------------------------------------------------
- REM SunAlign.BAS (Generic BASIC version)
- REM Calculates position of sun in sky, as azimuth (compass bearing
- REM measured clockwise from True North) and angle of elevation, as
- REM seen from any place on earth, on any date and any time.
- REM Also calculates alignment of a heliostat mirror.
- REM David Williams
- REM P.O. Box 48512
- REM 3605 Lakeshore Blvd. West
- REM Toronto, Ontario. M8W 4Y6
- REM Canada
- REM Original date 2007 Jul 07. This version 2007 Oct 07
- REM Note: For brevity, no error checks on user inputs
- CLS
- PRINT "Use negative numbers for opposite directions."
- INPUT "Observer's latitude (degrees North)"; LT
- INPUT "Observer's longitude (degrees East)"; LG
- INPUT "Date (M#,D#)"; Mth, Day
- INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
- INPUT "Time Zone (+
- hours from GMT/UT)"; TZN
PY = 4 * ATN(1): REM "PI" not assignable in some BASICs
DR = 180 / PY: REM degree/radian factor
W = 2 * PY / 365: REM earth's mean orbital speed in radians/day
C = -23.45 / DR: REM reverse angle of axial tilt in radians
ST = SIN(C): REM sine of reverse tilt
CT = COS(C): REM cosine of reverse tilt
E2 = 2 * .0167: REM twice earth's orbital eccentricity
SP = 12 * W: REM 12 days from December solstice to perihelion
D = INT(30.6 * ((Mth + 9) MOD 12) + 58.5 + Day) MOD 365
A = W * (D + 10): REM Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
C = (A - ATN(TAN(B) / CT)) / PY
ET = 720 * (C - INT(C + .5)): REM equation of time
REM in 720 minutes, earth rotates PI radians relative to sun
C = ST * COS(B)
EL = ATN(C / SQR(1 - C * C)) * DR: REM solar declination
AZ = 15 * (HR - TZN) + (MIN + ET) / 4 + LG: REM longitude diff
GOSUB 800
R = SQR(Y * Y + Z * Z)
AX = Y: AY = Z: GOSUB 710
A = AA + (90 - LT) / DR
Y = R * COS(A)
Z = R * SIN(A)
GOSUB 740
PRINT : REM AZ & EL are now sun's azimuth & elevation in degrees
IF EL < 0 THEN PRINT "Sun Below Horizon": END
R = AZ: GOSUB 870: PRINT "Sun's azimuth: "; R; " degrees"
R = EL: GOSUB 870: PRINT "Sun's elevation: "; R; " degrees"
PRINT
INPUT "Calculate heliostat mirror alignment (y/n)"; K$
IF K$ = "N" OR K$ = "n" THEN END
SX = X: SY = Y: SZ = Z
PRINT
INPUT "Azimuth of target direction (degrees)"; AZ
INPUT "Elevation of target direction (degrees)"; EL
GOSUB 800
X = X + SX: Y = Y + SY: Z = Z + SZ: GOSUB 740
PRINT : REM AZ & EL are now aim azimuth & elevation in degrees
PRINT "Mirror aim direction (perpendicular to surface):"
R = AZ: GOSUB 870: PRINT "Azimuth: "; R; " degrees"
R = EL: GOSUB 870: PRINT "Elevation: "; R; " degrees"
END
IF AX = 0 THEN AA = SGN(AY) * PY / 2: RETURN
AA = ATN(AY / AX): IF AX < 0 THEN AA = AA + PY
RETURN
AX = SQR(X * X + Y * Y): AY = Z: GOSUB 710
EL = AA * DR
AX = Y: AY = X: GOSUB 710
AZ = AA * DR
IF AZ < 180 THEN AZ = AZ + 180 ELSE AZ = AZ - 180
RETURN
E = EL / DR
A = AZ / DR
Z = SIN(E)
C = 0 - COS(E): REM Won't work without "0" in Liberty Basic
X = C * SIN(A)
Y = C * COS(A)
RETURN
R = INT(10 * R + .5): IF R = 3600 THEN R = 0
R = R / 10
RETURN
--------------------------------------------------
