Connect with us

heliostat program

Discussion in 'Home Power and Microgeneration' started by David Williams, Jul 10, 2007.

Scroll to continue with content
  1. 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

    ------------------------------------------------
     
Ask a Question
Want to reply to this thread or ask your own question?
You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.
Electronics Point Logo
Continue to site
Quote of the day

-