D
David Williams
- Jan 1, 1970
- 0
I recently wrote a QBasic program that calculates the position of the
sun in the sky, as azimuth and elevation, as seen from anywhere on
earth, at any time and date. It also calculates the orientation in which
a mirror should be pointed in order to reflect sunlight in any desired
direction, as is done by a heliostat. With a few hardware-specific
additions to make the computer control the mirror, this program could
be used to run a heliostat. In fact, it is very similar to the program
that I used to run a heliostat, years ago.
For anyone who may be interested, I'll append the program to this
message.
Incidentally, soon after I wrote it, I posted a version of this program
in another newsgroup. There was nothing wrong with that. It worked just
fine. But I have found a few ways to clean up and simplify the coding,
which are incorporated into this version.
Enjoy!
dow
--------------------------------------------------
' SunAlign.BAS
' Calculates position of sun in sky, as azimuth (true
' compass bearing) 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
' Toronto, Canada
' M8W 4Y6
' Initially dated 2007 Jul 07
' This version 2007 Jul 09
DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff$ (X)
DECLARE SUB D2P (X, Y, Z, AZ, EL)
DECLARE SUB P2D (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)
CONST PI = 3.1415926536#
CONST DR = 180 / PI ' degree / radian factor
CLS
PRINT "Use negative numbers for directions opposite to those shown."
INPUT "Observer's latitude (degrees North)"; LT
INPUT "Observer's longitude (degrees West)"; LG
INPUT "Date (M#,D#)"; Mth%, Day%
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
INPUT "Time Zone (+/- hours from GMT/UT)"; TZN
D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination
LD = 15 * (HR - TZN) + (MIN + ET) / 4 - LG 'longitude difference
CALL P2D(LD, DC, sX, sY, sZ)
RR = SQR(sY * sY + sZ * sZ)
CL = (90 - LT) / DR ' colatitude in radians
AA = Ang(sY, sZ) + CL
sY = RR * COS(AA)
sZ = RR * SIN(AA)
CALL D2P(sX, sY, sZ, sAZ, sEL)
PRINT
IF sEL < 0 THEN PRINT "Sun Below Horizon": END
PRINT "Sun's azimuth: "; ROff$(sAZ); " degrees"
PRINT "Sun's elevation: "; ROff$(sEL); " degrees"
PRINT
PRINT "Calculate heliostat mirror alignment? (y/n) ";
DO
K$ = UCASE$(INKEY$)
LOOP UNTIL K$ = "Y" OR K$ = "N"
PRINT K$
IF K$ = "Y" THEN
PRINT
INPUT "Azimuth of target direction (degrees)"; tAZ
INPUT "Elevation of target direction (degrees)"; tEL
CALL P2D(tAZ, tEL, tX, tY, tZ)
CALL D2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)
PRINT
PRINT "Mirror aim direction (perpendicular to surface):"
PRINT "Azimuth: "; ROff$(mAZ); " degrees"
PRINT "Elevation: "; ROff$(mEL); " degrees"
END IF
END
FUNCTION Ang (X, Y)
' calculates angle between positive X axis and vector to (X,Y)
IF X = 0 THEN
A = SGN(Y) * PI / 2
ELSE
A = ATN(Y / X)
IF X < 0 THEN A = A + PI
END IF
Ang = A
END FUNCTION
SUB D2P (X, Y, Z, AZ, EL)
' convert from X,Y,Z to AZ, EL
' Outputs in degrees
EL = Ang(SQR(X * X + Y * Y), Z) * DR
A = Ang(Y, X) * DR
IF A < 180 THEN AZ = A + 180 ELSE AZ = A - 180
END SUB
FUNCTION ET.Dec (D, F%) STATIC
' Calculates equation of time, in minutes, or solar declination,
' in degrees, on day number D of year. (D = 0 on January 1.)
' F% selects function: True (non-zero) for Equation of Time,
' False (zero) for Declination.
' STATIC means variables are preserved between calls of function
' This version assumes PI and DR (180/PI) are already initialized
IF W = 0 THEN ' first call, initialize constants
W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
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
SP = 12 * W ' 12 days from December solstice to perihelion
D1 = -1 ' holds last D. Saves time if D repeated for both functions
END IF
IF D <> D1 THEN ' new value of D
A = W * (D + 10) ' Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
D1 = D
END IF
IF F% THEN ' equation of time calculation
C = (A - ATN(TAN(B) / CT)) / PI
ET.Dec = 720 * (C - INT(C + .5))
' in 720 minutes, earth rotates PI radians relative to sun
ELSE ' declination calculation
C = ST * COS(B)
ET.Dec = ATN(C / SQR(1 - C * C)) * DR
' arcsine of C in degrees. ASN not directly available in QBasic
END IF
END FUNCTION
SUB P2D (AZ, EL, X, Y, Z)
' convert from AZ,EL to X,Y,Z
' Inputs in degrees
E = EL / DR
A = AZ / DR
Z = SIN(E)
C = -COS(E)
X = C * SIN(A)
Y = C * COS(A)
END SUB
FUNCTION ROff$ (X)
' neatly rounds number to one place of decimals
S$ = LTRIM$(STR$(INT(10 * ABS(X) + .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$
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1)
END FUNCTION
------------------------------------------------
sun in the sky, as azimuth and elevation, as seen from anywhere on
earth, at any time and date. It also calculates the orientation in which
a mirror should be pointed in order to reflect sunlight in any desired
direction, as is done by a heliostat. With a few hardware-specific
additions to make the computer control the mirror, this program could
be used to run a heliostat. In fact, it is very similar to the program
that I used to run a heliostat, years ago.
For anyone who may be interested, I'll append the program to this
message.
Incidentally, soon after I wrote it, I posted a version of this program
in another newsgroup. There was nothing wrong with that. It worked just
fine. But I have found a few ways to clean up and simplify the coding,
which are incorporated into this version.
Enjoy!
dow
--------------------------------------------------
' SunAlign.BAS
' Calculates position of sun in sky, as azimuth (true
' compass bearing) 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
' Toronto, Canada
' M8W 4Y6
' Initially dated 2007 Jul 07
' This version 2007 Jul 09
DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff$ (X)
DECLARE SUB D2P (X, Y, Z, AZ, EL)
DECLARE SUB P2D (AZ, EL, X, Y, Z)
DECLARE FUNCTION Ang (X, Y)
CONST PI = 3.1415926536#
CONST DR = 180 / PI ' degree / radian factor
CLS
PRINT "Use negative numbers for directions opposite to those shown."
INPUT "Observer's latitude (degrees North)"; LT
INPUT "Observer's longitude (degrees West)"; LG
INPUT "Date (M#,D#)"; Mth%, Day%
INPUT "Time (HH,MM) (24-hr format)"; HR, MIN
INPUT "Time Zone (+/- hours from GMT/UT)"; TZN
D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination
LD = 15 * (HR - TZN) + (MIN + ET) / 4 - LG 'longitude difference
CALL P2D(LD, DC, sX, sY, sZ)
RR = SQR(sY * sY + sZ * sZ)
CL = (90 - LT) / DR ' colatitude in radians
AA = Ang(sY, sZ) + CL
sY = RR * COS(AA)
sZ = RR * SIN(AA)
CALL D2P(sX, sY, sZ, sAZ, sEL)
IF sEL < 0 THEN PRINT "Sun Below Horizon": END
PRINT "Sun's azimuth: "; ROff$(sAZ); " degrees"
PRINT "Sun's elevation: "; ROff$(sEL); " degrees"
PRINT "Calculate heliostat mirror alignment? (y/n) ";
DO
K$ = UCASE$(INKEY$)
LOOP UNTIL K$ = "Y" OR K$ = "N"
PRINT K$
IF K$ = "Y" THEN
INPUT "Azimuth of target direction (degrees)"; tAZ
INPUT "Elevation of target direction (degrees)"; tEL
CALL P2D(tAZ, tEL, tX, tY, tZ)
CALL D2P(sX + tX, sY + tY, sZ + tZ, mAZ, mEL)
PRINT "Mirror aim direction (perpendicular to surface):"
PRINT "Azimuth: "; ROff$(mAZ); " degrees"
PRINT "Elevation: "; ROff$(mEL); " degrees"
END IF
END
FUNCTION Ang (X, Y)
' calculates angle between positive X axis and vector to (X,Y)
IF X = 0 THEN
A = SGN(Y) * PI / 2
ELSE
A = ATN(Y / X)
IF X < 0 THEN A = A + PI
END IF
Ang = A
END FUNCTION
SUB D2P (X, Y, Z, AZ, EL)
' convert from X,Y,Z to AZ, EL
' Outputs in degrees
EL = Ang(SQR(X * X + Y * Y), Z) * DR
A = Ang(Y, X) * DR
IF A < 180 THEN AZ = A + 180 ELSE AZ = A - 180
END SUB
FUNCTION ET.Dec (D, F%) STATIC
' Calculates equation of time, in minutes, or solar declination,
' in degrees, on day number D of year. (D = 0 on January 1.)
' F% selects function: True (non-zero) for Equation of Time,
' False (zero) for Declination.
' STATIC means variables are preserved between calls of function
' This version assumes PI and DR (180/PI) are already initialized
IF W = 0 THEN ' first call, initialize constants
W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
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
SP = 12 * W ' 12 days from December solstice to perihelion
D1 = -1 ' holds last D. Saves time if D repeated for both functions
END IF
IF D <> D1 THEN ' new value of D
A = W * (D + 10) ' Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
D1 = D
END IF
IF F% THEN ' equation of time calculation
C = (A - ATN(TAN(B) / CT)) / PI
ET.Dec = 720 * (C - INT(C + .5))
' in 720 minutes, earth rotates PI radians relative to sun
ELSE ' declination calculation
C = ST * COS(B)
ET.Dec = ATN(C / SQR(1 - C * C)) * DR
' arcsine of C in degrees. ASN not directly available in QBasic
END IF
END FUNCTION
SUB P2D (AZ, EL, X, Y, Z)
' convert from AZ,EL to X,Y,Z
' Inputs in degrees
E = EL / DR
A = AZ / DR
Z = SIN(E)
C = -COS(E)
X = C * SIN(A)
Y = C * COS(A)
END SUB
FUNCTION ROff$ (X)
' neatly rounds number to one place of decimals
S$ = LTRIM$(STR$(INT(10 * ABS(X) + .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$
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1)
END FUNCTION
------------------------------------------------
