Space Forum / Space Science / January 2008

Well... Almost astronomy, anyway.
 
In an effort to get this newsgroup back on topic, and also to show how
simple astronomical mathematics can be, here's a QBasic computer
program that calculates where the sun is in the sky, as seen from
anywhere on the earth, at any time on any date. It also calculates how
a mirror should be aligned if it is to reflect sunlight in a desired
direction. This is what is done by "heliostat" mirrors, which are used
in various kinds of solar-energy applications, for example.
 
The calculations are a bit approximate. They assume that the earth is a
perfect sphere, for example, and they ignore leap-years. Even so, the
results are accurate to a smallish fraction of a degree, which is
plenty good enough for solar-energy applications.
 
I don't think I have posted this in this newsgroup before, but people
who follow other 'groups may have seen it.
 
The thing I like the most is how elegant and simple the calculations
turn out to be. There are only about twenty lines of mathematical code
in the whole thing.
 
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 2007 Dec 23
 
' 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 FUNCTION DegOut$ (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
 
  PRINT "Sun's azimuth: "; DegOut$(sAZ)
  PRINT "Sun's elevation: "; DegOut$(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):"
  PRINT "Azimuth: "; DegOut$(mAZ)
  PRINT "Elevation: "; DegOut$(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
 
FUNCTION DegOut$ (X)
' converts angle to degrees rounded to one place of decimals
 
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$
DegOut$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1) + " degrees"
 
END FUNCTION
 
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
 
--------------------------------------------------------