Re: heliostat program
- From: "JJ" <stackoverflow@xxxxxxxxxxx>
- Date: Thu, 12 Jul 2007 03:01:32 +0100
"David Williams" <david.williams@xxxxxxxxxx> wrote in message news:1184201508.877.0@xxxxxxxxxxxxx
A few days ago, I wrote a QBasic program that calculates the positionLooking over this, I note that there are a couple of calls to
of the sun in the sky, as azimuth (true compass bearing) and angle of
elevation, as seen from anywhere on earth, at any time on any date. It
also calculates the required alignment of a mirror if it is to reflect
sunlight in any desired direction, as is done in a heliostat. With some
hardware-specific additions, to make the computer control motors that
turn the mirror, this program could be the basis of the software for a
computer-controlled heliostat. In fact, it is similar to a VIC 20
program of mine that I used for that purpose, years ago.
The calculations contain some minor approximations, but they should be
accurate enough for all purposes except precise astronomy.
In case anyone would like to have it, here it is.
dow
---------------------------------------------
' SunAlign
' 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
subroutines where some of the variables being passed seem
to not been initialised.
Does QBasic assign zero to such variables in these circumstances?
e.g
CALL P2D(LD, DC, sX, sY, sZ)
CALL D2P(sX, sY, sZ, sAZ, sEL)
JJ
.
- Follow-Ups:
- Re: heliostat program
- From: Stephen J. Rush
- Re: heliostat program
- References:
- heliostat program
- From: David Williams
- heliostat program
- Prev by Date: heliostat program
- Next by Date: Re: heliostat program
- Previous by thread: heliostat program
- Next by thread: Re: heliostat program
- Index(es):
Relevant Pages
|
