' 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