' 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 Sep 02
' 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)
AZ = Ang(Y, X)
IF AZ < 0 THEN AZ = AZ + 2 * 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