Using solar energy frequently involves calculations. This page is to show how various calculations can be done. If you have calculated anything to do with solar energy and are prepared to share your method with other people, please post it here. Just click on the "edit" button, above, and add your material near the bottom of the page, just above the "External Links" line. Put a couple of "equals" signs (=) before and after the title of your addition. This will make Appropedia format it properly. If the title is done correctly, your item will automatically be added to the table of contents at the top of the page. If you have any problems with formatting and the like, don't worry. Someone will tidy it up for you. Or, if you're uncomfortable with editing the article, go to the "discussion" page and put your calculation there. Someone will move it into the main article for you.
Code to calculate Equation of Time and Solar Declination[edit | edit source]
The Equation of Time is the difference between the time that is shown by a sundial and the time that is shown by a clock that is set to local mean time. Essentially, the Equation of Time shows the longitude of the sun in the sky, relative to its mean position. The Solar Declination is the sun's latitude. Both of these quantities must be known for many purposes related to solar energy.
Computer program[edit | edit source]
The following subroutine, here written in QBasic but easily translatable to other languages, calculates the Equation of Time and the Solar Declination on any day of the year. It is quite accurate. Its Root-Mean-Square error for the Equation of Time is only 3.7 seconds. The Declination is calculated with errors that are always small compared with the angular radius of the sun as seen from the earth (about 0.25 degrees).
A program that includes the routine is ETIMSDEC, which is listed below. Anyone who wants further information should look at the Sun related texts page. The article titled "The Latitude and Longitude of the Sun", which I wrote several years ago, describes the astronomical logic behind the routine.
DOwenWilliams 10:19, 18 May 2011 (PDT) David Williams
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 IF PI = 0 THEN ' first call, initialize constants PI = 4 * ATN(1) W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day DR = 180 / PI ' degree/radian factor 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)) ' this is value of equation of time ' 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 ' this is value of declination ' arcsine of C in degrees. ASN not directly available in QBasic END IF END FUNCTION
Alternative calculation[edit | edit source]
The following calculation is essentially the same as the one above, with minor changes. Angles are in degrees instead of radians. Some explanation follows each line. And, although it can be written as a computer program, the calculation can be done with a pocket calculator or similar tool.
Angles are in degrees. Asterisks mean multiplication. Slashes mean division. The conventional order of operations applies.
W is the Earth's mean angular orbital velocity in degrees per day.
D is the date, in days starting at zero on January 1 (i.e. the days part of the ordinal date -1). 10 is the approximate number of days from the December solstice to January 1. "A" is the angle the earth would move on its orbit at its average speed from the December solstice to date D.
B is the angle the Earth moves from the solstice to date D, including a first-order correction for the Earth's orbital eccentricity, 0.0167. The number 2 is the number of days from January 1 to the date of the Earth's perihelion. This expression for B can be simplified by combining constants to: .
C is the difference between the angles moved at mean speed, and at the corrected speed projected onto the equatorial plane, and divided by 180 to get the difference in "half turns". The number 23.44 is the obliquity (tilt) of the Earth's axis in degrees. The subtraction gives the conventional sign to the equation of time. For any given value of x, arctan(x) (sometimes written as tan−1x) has multiple values, differing from each other by integer numbers of half turns. The value generated by a calculator or computer may not be the appropriate one for this calculation. This may cause C to be wrong by an integer number of half turns. The excess half turns are removed in the next step of the calculation:
EoT is the equation of time in minutes. The expression nint(C) means the nearest integer to C. On a computer, it can be programmed, for example, as INT(C+0.5). It is 0, 1, or 2 at different times of the year. Subtracting it leaves a small positive or negative fractional number of half turns, which is multiplied by 720, the number of minutes (12 hours) that the Earth takes to rotate one half turn relative to the Sun, to get the equation of time.
Compared with published values, this calculation has a Root Mean Square error of only 3.7 seconds of time. The greatest error is 6.0 seconds.
Addendum about solar declination[edit | edit source]
The value of B in the above calculation is an accurate value for the Sun's ecliptic longitude (shifted by 90 degrees), so the solar declination becomes readily available:
which is accurate to within a fraction of a degree.
Program to calculate sun's position and heliostat mirror alignment[edit | edit source]
This program, written in QBasic but easily translatable to other programming languages, calculates the position of the sun in the sky, as angles of azimuth (true compass bearing) and elevation, as seen from any place on the earth, at any time on any date. It also calculates the correct alignment of a heliostat mirror to reflect sunlight onto any specified target. It is accurate to within 0.25 degree.
For detailed explanations of how the program works, see the Sun related texts page of this wiki.
I am the author of this program. My name and address are on it so anyone who has questions or problems can contact me. The program is in the public domain, not copyright. Anyone is free to use it, or any part of it, for any purpose, commercial or non-commercial.
DOwenWilliams 07:46, 27 April 2012 (PDT)
' SunAlign.BAS (Version for QBasic and similar dialects) ' Calculates position of sun in sky, as azimuth (compass bearing) ' and angle of elevation, and orientation of heliostat mirror. ' David Williams ' P.O. Box 48512 ' Long Branch P.O. ' Toronto, Ontario. M8W 1P8 ' Canada ' Initially dated 2007 Jul 07 ' This version 2008 Nov 03 ' 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 PY2 = PY + PY CONST DR = 180 / PY ' degree / radian factor W = PY2 / 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 PRINT "SunAlign" PRINT PRINT "Calculates position of sun in sky, as angles of azimuth and" PRINT "elevation, as seen from any place on earth, on any date and at" PRINT "any time. Also calculates the alignment of a heliostat mirror." PRINT "Could be modified to make software to drive a heliostat or" PRINT "sun tracker." PRINT PRINT "Azimuths in degrees (0 to 359.9) clockwise from True North." PRINT "0 = North, 90 = East, 180 = South, 270 = West." PRINT PRINT "Angles of elevation in degrees upward from horizontal." PRINT "Negative elevations (downward) allowed for target and mirror." PRINT Menu: PRINT 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 east-west 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" GOTO Menu FUNCTION Ang (X, Y) ' calculates angle at origin 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 + PY2 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 * VAL(S$) < 0 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
Estimating sunrise and sunset data, using analemma[edit | edit source]
If it is marked with the dates on which the Sun is at various points on it, the analemma summarises the apparent movements of the Sun, relative to its mean position, throughout the year. A diagram of the analemma, with equal scales in both north-south and east-west directions, can be used as a tool to estimate quantities such as the times of sunrise and sunset, which depend on the Sun's position. Generally, making these estimates depends on visualizing the analemma as a rigid structure in the sky, which moves around the Earth at constant speed so it rises and sets once a day, with the Sun slowly moving around it once a year. Some approximations are involved in the process, chiefly the use of a plane diagram to represent things on the celestial sphere, and the use of drawing and measurement instead of numerical calculation. Because of these, the estimates are not perfectly precise, but they are usually good enough for practical purposes. Also, they have instructional value, showing in a simple visual way how the times of sunrises and sunsets vary.
Earliest and latest sunrise and sunset[edit | edit source]
The analemma can be used to find the dates of the earlest and latest sunrises and sunsets of the year. These do not occur on the dates of the solstices.
If we visualize the analemma in the eastern sky, with its lowest point just above the horizon, that point has just risen above the horizon. If the Sun were at that point, sunrise would have just occurred. This would be the latest sunrise of the year, since all other points on the analemma would rise earlier. Therefore the date when the Sun is at this lowest point is the date of the latest sunrise. Similarly, when the Sun is at the highest point on the analemma, the earliest sunrise of the year will occur. Likewise, at sunset, the earliest sunset will occur when the Sun is at its lowest point on the analemma when it is close to the western horizon, and the latest sunset when it is at the highest point.
None of these points is exactly at one of the ends of the analemma, where the Sun is at a solstice. As seen from north-mid-temperate latitudes, the earliest sunset occurs some time before the December solstice - typically a week or two before it - and the latest sunrise happens a week or two after the solstice. Thus, the darkest evening occurs in early to mid-December, but the mornings keep getting darker until about the New Year.
The exact dates are those on which the Sun is at the points where the horizon is tangential to the analemma, which in turn depend on how much the analemma, or the north-south meridian passing through it, is tilted from the vertical. This angle of tilt is essentially the co-latitude (90 degrees minus the latitude) of the observer. Calculating these dates numerically is complex, but they can be estimated fairly accurately by placing a straight-edge, tilted at the appropriate angle, tangential to a diagram of the analemma, and reading the dates when the Sun is at the positions of contact. In temperate latitudes, the dates get further from the solstices as the latitude decreases. In near-equatorial latitudes, the situation is more complex. The analemma lies almost horizontal, and there are two widely separated dates in the year when the Sun rises earlier than on adjoining dates, and so on.
Times of sunrise and sunset[edit | edit source]
A similar geometrical method, based on the analemma, can be used to find the times of sunrise and sunset at any place on the Earth (except within or close to the Arctic or Antarctic Circle), on any date.
The "origin" of the analemma, where the solar declination and the equation of time are both zero, rises and sets at 6 am and 6 pm local mean time on every day of the year, irrespective of the latitude of the observer. (This estimation does not take account of atmospheric refraction.) If the analemma is drawn in a diagram, tilted at the appropriate angle for an observer's latitude (as described above), and if a horizontal line is drawn to pass through the position of the Sun on the analemma on any given date, then at sunrise this line represents the horizon. The origin moves along the celestial equator at a speed of 15 degrees per hour, the speed of the Earth's rotation. The distance along the equator from the point where it intersects the horizon to the position of the origin of the analemma at sunrise is the distance the origin moves between 6 am and the time of sunrise on the given date. Measuring the length of this equatorial segment therefore gives the difference between 6 am and the time of sunrise.
The measurement should, of course, be done on the diagram, but it should be expressed in terms of the angle that would be subtended at an observer on the ground by the corresponding distance in the analemma in the sky. It can be useful to compare it with the length of the analemma, which subtends 47 degrees. Thus, for example, if the length of the equatorial segment on the diagram is 0.2 times the length of the analemma on the diagram, then the segment in the celestial analemma would subtend 0.2 x 47 or 9.4 degrees at the observer on the ground. The angle, in degrees, should be divided by 15 to get the time difference in hours between sunrise and 6 am. The sign of the difference is clear from the diagram. If the horizon line at sunrise passes above the origin of the analemma, the Sun rises before 6 am, and vice versa.
The same technique can be used, mutatis mutandis, to estimate the time of sunset. Note that the estimated times are in local mean time. Corrections must be applied to convert them to standard time or daylight saving time. These corrections will include a term that involves the observer's longitude, so both his latitude and longitude affect the final result.
Azimuths of sunrise and sunset[edit | edit source]
The azimuths (true compass bearings) of the points on the horizon where the Sun rises and sets can be easily estimated, using the same diagram as is used to find the times of sunrise and sunset, as described above.
The origin of the analemma rises due east, and sets due west. Therefore, on the diagram, the point where the horizon intersects the equator, which is where the origin crosses the horizon, must represent due east or west. The point where the Sun is at sunrise or sunset represents the direction of sunrise or sunset. Simply measuring the distance along the horizon between these points, in angular terms, gives the angle between due east or west and the direction of sunrise or sunset.
Equation of Time, Declination, and Analemma Program[edit | edit source]
The following QBasic program displays graphs of the Equation of Time and the Solar Declination, showing comparisons with published data. It also draws the analemma, with different colours used to show the parts of the analemma where the sun is in different months.
' ETIMSDEC.BAS
' David Williams, 2003
' P.O. Box 48512
' 3701 Lakeshore Blvd. West
' Toronto, Ontario. M8W 1P8
' Canada
' This version dated 2007 Jun 13
DECLARE FUNCTION ET.Dec (D, F%)
DECLARE FUNCTION ROff$ (X)
DECLARE SUB Waitkey (X%)
CL0 = 7: CL1 = 12: CL2 = 15 ' colours used
SCREEN 12
COLOR CL2
CLS
PRINT "ETIMSDEC.BAS"
PRINT
PRINT "Shows Equation of Time and Solar Declination calculations"
PRINT "(performed in function ET.Dec) and compares their results"
PRINT "graphically with published values, showing close agreement."
PRINT
PRINT "These calculations can be used in the programs of computer-"
PRINT "controlled solar-energy equipment, such as sun trackers and"
PRINT "heliostats."
PRINT
PRINT "Equation of Time is difference between Solar Time and Mean"
PRINT "Time. Sundials show solar time. Clocks show mean time. The"
PRINT "Equation of Time is related to solar longitude, relative to"
PRINT "its mean value (the value it would have if the sun's apparent"
PRINT "motion were at uniform speed). One degree westward of solar"
PRINT "longitude corresponds to four minutes in positive direction in"
PRINT "the Equation of Time."
PRINT
PRINT "Solar Declination is latitude of sun in celestial coordinates."
PRINT
' initialize array used in graphing routine
DIM ML(1 TO 13)
FOR X = 1 TO 13
READ ML(X)
NEXT
DATA 33,29,33,31,33,31,33,33,31,33,31,33,0
' Month lengths adjusted for numbers of pixels on screen
DO
COLOR CL0
PRINT "1. Find Equation of Time and Solar Declination on given date"
PRINT "2. Draw Equation of Time Graph"
PRINT "3. Draw Solar Declination Graph"
PRINT "4. Draw Solar Analemma"
PRINT "5. Quit program"
PRINT
PRINT "Which (1 - 5)? ";
WHILE INKEY$ <> "": WEND
DO
K$ = INKEY$
LOOP UNTIL K$ >= "1" AND K$ <= "5"
PRINT K$
SELECT CASE VAL(K$)
CASE 1: GOSUB Calc
CASE 2, 3: GOSUB Graph
CASE 4: GOSUB Analemma
CASE 5: EXIT DO
END SELECT
LOOP
GOTO WayOut
Calc:
PRINT
COLOR CL1
DO
INPUT "Date (M#,D#)"; Mth%, Day%
LOOP UNTIL Mth% > 0 AND Mth% < 13 AND Day% > 0 AND Day% < 32
D = INT(30.6 * ((Mth% + 9) MOD 12) + 58.5 + Day%) MOD 365
' day of 365-day year. Jan 1 = 0. Dec 31 = 364
PRINT
PRINT "Day number"; D + 1; "of year" ' conventional range, 1 - 365
ET = ET.Dec(D, 1) ' Equation of Time
DC = ET.Dec(D, 0) ' Declination
PRINT
PRINT "Equation of Time: "; ROff$(ET); " minutes"
PRINT "Solar Declination: "; ROff$(DC); " degrees"
PRINT
RETURN ' to menu
Graph:
F% = (K$ = "2") ' true if Equation of Time selected
' Set up graph
CLS
IF F% THEN
PRINT TAB(30); "EQUATION OF TIME"
LOCATE 3, 15
PRINT "Graph shows difference in minutes between clock"
PRINT TAB(15); "and sundial time. Positive difference means"
PRINT TAB(15); "sundial is ahead of clock, and vice versa."
LOCATE 24, 15
COLOR CL1
PRINT "Graph shows results of calculation."
COLOR CL2
PRINT TAB(15); "Crosses show published values of Equation of Time."
RESTORE ETData
ELSE
PRINT TAB(25); "Declination of Sun, in Degrees"
LOCATE 3, 10
COLOR CL1
PRINT "Graph shows calculated function. ";
COLOR CL2
PRINT "Crosses show published values."
RESTORE DeclnData
END IF
READ ML%, VT%, VB%, LL%, UL%
COLOR CL0
LOCATE 16, 67
PRINT "-="
FOR T = LL% TO UL% STEP 5
LINE (137, 247 - 6.4 * T)-(530, 247 - 6.4 * T)
LOCATE 16 - T / 2.5, 14
PRINT RIGHT$(" " + STR$(T), 3);
IF T = 0 THEN PRINT " ="
NEXT
LOCATE ML%, 20
PRINT "Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec"
' Do graphing
T = 0
D = 0
MC = 0
FOR M = 0 TO 384
M1 = 147 + M
IF M = T THEN
' Draw vertical line
COLOR CL0
LINE (M1, VT%)-(M1, VB%)
MC = MC + 1
T = T + ML(MC)
MS = 0
ELSE
' Plot point(s)
COLOR CL1
YY = 247 - 6.4 * ET.Dec(D, F%) ' function used
IF POINT(M1, YY) <> CL2 THEN PSET (M1, YY) ' don't overwrite cross
IF MS = 0 OR MS = 10 OR MS = 20 THEN
' Plot cross representing published value
COLOR CL2
READ E ' published value
E = FIX(E) + (E - FIX(E)) * 5 / 3 ' see note after E of T data
YY = 247 - 6.4 * E
LINE (M1 - 2, YY)-(M1 + 2, YY) ' lines make
LINE (M1, YY - 2)-(M1, YY + 2) ' cross
END IF
MS = MS + 1
D = D + 1
END IF
IF ML(MC) = 33 AND M = T - 18 THEN M = M + 1
NEXT
Waitkey 25
RETURN ' to menu
ETData:
DATA 7, 97, 343, -15, 20: ' plotting parameters for E. of T. graph
DATA -3.12,-7.38,-11.08,-13.33,-14.19,-13.49,-12.34,-10.18,-7.28
DATA -4.08,-1.16,1,2.51,3.4,3.34,2.25,.39,-1.28
DATA -3.33,-5.16,-6.15,-6.16,-5.14,-3.16,-.12,3.08,6.4
DATA 10.05,13.02,15.12,16.2,16,14.16,11.11,7.02,2.13
' Above data show published values of Equation of Time on
' 1st, 11th and 21st days of months. Decimal point separates
' minutes, to left, from seconds, to right, e.g. 3.16 means
' 3 minutes and 16 seconds. Seconds part is multiplied by 5/3
' to get fractional minutes after READ instruction.
DeclnData:
DATA 5, 65, 407, -25, 25: ' plotting parameters for Decl'n graph
DATA -23.04,-21.56,-20.05,-17.2,-14.18,-10.52,-7.49,-3.57,0
DATA 4.18,8.07,11.39,14.54,17.43,20.04,21.58,23.02,23.26
DATA 23.09,22.11,20.36,18.1,15.27,12.19,8.3,4.48,.57
DATA -2.57,-6.48,-10.29,-14.14,-17.15,-19.47,-21.43,-22.57,-23.26
' Above data show published values (in degrees and minutes) of solar
' declination on 1st, 11th and 21st days of month. See note after
' Equation of Time data.
' Published values are from the book:
' Sundials, Their Theory and Construction
' by: Albert E. Waugh
' Dover Publications, New York, 1973
' ISBN 0-486-22947-5
Analemma:
CLS
RESTORE Anadata
PRINT "SOLAR ANALEMMA"
PRINT
PRINT
PRINT "Shows apparent path of sun,"
PRINT "through the year, relative"
PRINT "to its mean position"
PRINT "(marked by cross)."
LOCATE 3, 60
PRINT "N"
LOCATE 7, 60
PRINT "S"
LOCATE 5, 56
PRINT "E W"
LINE (475, 52)-(475, 90) ' compass lines
LINE (456, 71)-(494, 71)
LOCATE 9, 1
PRINT TAB(53); "Looking at sky,"
PRINT TAB(53); "if North is upward,"
PRINT TAB(53); "East is to left."
LOCATE 24, 1
PRINT "For dimensions, see graphs"
PRINT "of Solar Declination and"
PRINT "Equation of Time."
PRINT "+4 minutes in Equation of Time"
PRINT "equal 1 degree Westward."
LOCATE 15, 1
PRINT TAB(50); "Celestial"
PRINT TAB(50); "Equator"
X0 = 320 ' origin coordinates
Y0 = 240
LINE (X0, Y0 + 5)-(X0, Y0 - 5) ' cross at (0,0)
LINE (X0 + 5, Y0)-(X0 - 5, Y0)
FOR S = -1 TO 1 STEP 2
LINE (X0, Y0 + S * 210)-(X0, Y0 + S * 225) ' lines marking N-S
LINE (X0 + S * 30, Y0)-(X0 + S * 55, Y0) ' and E-W
NEXT
LOCATE 10, 1
PRINT "Month colours:";
M = 0 ' month number
DM = 0 ' day of month
LM = 0 ' length of month
' draw analemma
FOR DY = 0 TO 364 ' day of year
IF DM = LM THEN ' start new month
M = M + 1
C = M XOR 14 * (M AND 1) ' twiddle bits to get good contrasts
IF C = 8 THEN C = 14 ' yellow better than dull grey
COLOR C
READ LM, MN$
PRINT TAB(17); MN$ ' month name in correct colour
DM = 0 ' reset day of month
END IF
DM = DM + 1
FOR D = DY TO DY + .9 STEP .25 ' 4 points per day to avoid gaps
E = ET.Dec(D, 1) ' equation of time
L = ET.Dec(D, 0) ' solar declination (latitude)
PSET (X0 + 2.5 * E, Y0 - 10 * L) ' 4:1 ratio. 4 mins = 1 degree
NEXT
NEXT
Waitkey 48
RETURN ' to menu
Anadata:
DATA 31,Jan,28,Feb,31,Mar,30,Apr,31,May,30,Jun
DATA 31,Jul,31,Aug,30,Sep,31,Oct,30,Nov,31,Dec
' month lengths and names
WayOut:
COLOR CL2
CLS
PRINT "These calculations of Equation of Time and Solar Declination"
PRINT "are simplified and approximate. However, they are quite good."
PRINT "All differences between calculated and published values of"
PRINT "the solar declination are small compared with the angular"
PRINT "size of the sun in the sky. Similarly, differences between"
PRINT "calculated and published values of the equation of time are"
PRINT "small compared with the time the sun takes to traverse its"
PRINT "own diameter as it moves across the sky. The non-zero size"
PRINT "of the sun, rather than any inaccuracies of calculation,"
PRINT "is the limiting factor on how accurately the sun can be"
PRINT "tracked in solar energy applications, using this routine."
PRINT
PRINT "Leap years cause minor fluctuations in the Equation of Time"
PRINT "and Solar Declination on given dates of the year. However,"
PRINT "the errors that are introduced by ignoring these fluctuations"
PRINT "are always smaller than the apparent radius of the sun. This"
PRINT "program therefore ignores leap years. February 29, when it"
PRINT "occurs, is considered to be the same day as March 1."
PRINT
PRINT "The Equation of Time is treated here as the correction to be"
PRINT "subtracted from a sundial reading to get local mean (clock)"
PRINT "time. It is therefore positive when sundial is ahead of clock."
PRINT "This is the usual sign convention, but the opposite usage is"
PRINT "sometimes found. Take care when comparing values of the"
PRINT "Equation of Time from different sources."
Waitkey 1
PRINT "Note that atmospheric refraction of light affects the apparent"
PRINT "position of the sun in the sky when it is close to the"
PRINT "horizon. Sunlight is then too weak to be used for most"
PRINT "solar-energy applications, so refraction is not usually"
PRINT "included in calculations related to solar energy. This program"
PRINT "does not take account of refraction."
PRINT
PRINT "Comments in the function ET.Dec show how quantities in it are"
PRINT "related to astronomical parameters. Over long periods of time"
PRINT "(many decades or more), some of these parameters will change"
PRINT "significantly. The function should be modified accordingly if"
PRINT "it is to be used over very long periods of time."
Waitkey 1
END
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
IF PI = 0 THEN ' first call, initialize constants
PI = 4 * ATN(1)
W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
DR = 180 / PI ' degree/radian factor
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
FUNCTION ROff$ (X)
' neatly rounds number to one place of decimals
S$ = LTRIM$(STR$(INT(10 * ABS(X) + .5)))
IF LEN(S$) = 1 THEN S$ = "0" + S$
IF VAL(S$) THEN S$ = MID$("-x+", SGN(X) + 2, 1) + S$
ROff$ = LEFT$(S$, LEN(S$) - 1) + "." + RIGHT$(S$, 1)
END FUNCTION
SUB Waitkey (X%)
' prints prompt at position X% on bottom row, then waits
COLOR 15
LOCATE 30, X%
PRINT "*** Press any key to continue ***";
WHILE INKEY$ <> "": WEND
WHILE INKEY$ = "": WEND
CLS
END SUB
Formulae about paraboloidal dishes[edit | edit source]
Dimensions of a dish[edit | edit source]
The dimensions of a symmetrical paraboloidal dish are related by the equation: where is the focal length, is the depth of the dish (measured along the axis of symmetry from the vertex to the plane of the rim), and is the radius of the rim. Of course, they must all be in the same units. If two of these three quantities are known, this equation can be used to calculate the third.
A more complex calculation is needed to find the diameter of the dish measured along its surface. This is sometimes called the "linear diameter". It is the diameter of a flat, circular sheet of material, usually metal, which is the right size to be cut and bent to make the dish. Two intermediate results are useful in the calculation: (or equivalently and where and are defined as above. The linear diameter, is then given by: where means the natural logarithm of , i.e. its logarithm to base "e". (This calculation is exactly accurate. Computer programs that approximate a parabola with a short sequence of straight line segments are less accurate. This includes some published programs.)
The volume of the dish, the amount of liquid it could hold if the rim were horizontal and the vertex at the bottom, is given by or equivalently or where the symbols are defined as above. This last version can be compared with the well-known formulae for the volumes of a cylinder and a cone. Of course, is the aperture area of the dish, the area enclosed by the rim, which is proportional to the amount of sunlight the reflector dish can intercept.
The surface area of a paraboloidal dish can be found using the area formula for a surface of revolution, which gives , providing .
Paraboloids made by rotating liquids[edit | edit source]
Paraboloids can be made quite simply by rotating open-topped containers which hold liquids. The top surface of a liquid which is being rotated at constant speed around a vertical axis naturally takes the form of a paraboloid.
The focal length of the paraboloid is related to the angular speed at which the liquid is rotated by the equation: 2f=g, where f is the focal length, is the rotation speed, and g is the acceleration due to gravity. They must be in compatible units so, for example, f can be in metres, in radians per second, and g in metres per second-squared. The angle unit in must be radians. 1 radian per second is about 9.55 rotations per minute (RPM). On the Earth's surface, g is about 9.81 metres per second-squared. Putting these numbers into the equation produces the approximation: , where is the focal length in metres, and s is the rotation speed in RPM.
Focus-balanced reflectors[edit | edit source]
Paraboloidal reflectors that have their centres of mass coincident with their focal points are useful. If the paraboloidal reflector is axially symmetrical and is made of material of uniform thickness, its centre of mass coincides with its focus if the depth of the reflector, measured along its axis of symmetry from the vertex to the plane of the rim, is 1.8478 times its focal length. The radius of the rim of the reflector is 2.7187 times the focal length. The angular radius of the rim, as seen from the focal point, is 72.68 degrees. For details of how these numbers can be calculated, see the "Main article".
Iterative calculations[edit | edit source]
There are situations in which a calculation is impracticably difficult or impossible to do, even by a professional mathematician. Often, in these situations, a numerical answer can be obtained by using a computer program that performs a large number of iterations around a loop. Using these techniques requires the availability of a computer programming language (it makes little difference which one), and some proficiency in using it. Obtaining these is very useful.
Numerical integration[edit | edit source]
Integration is an operation that can be extremely difficult or impossible. However, integration is really just summation of an infinite number of infinitessimal terms. By making the terms larger than infinitessimal, but still very small, a computer program can make the machine add a huge number of terms to come up with an acceptably accurate value for the definite integral. An example is shown in Focus-balanced paraboloid. The integration was impossibly difficult, at least for the writer. Writing the computer program was quite easy, and the answer was accurate to at least ten significant digits.
Binary searches to invert functions[edit | edit source]
In "Dimensions of a dish", above, an expression is given that allows the "linear diameter" of a paraboloidal dish, which is represented by L, to be calculated from its focal length, F, and the radius of its rim, R. Doing the calculation in that direction is easy, using a pocket calculator. However, inverting it to calculate F from L and R, for example, is extremely difficult. One approach is to program a computer to perform a binary or bisection method search for the correct value of F. Using the given R and a trial value of F which is at the middle of an assumed range of possibility, the machine calculates a value for L. Depending on whether the calculated L is too high or too low, compared with the given one, the trial F is assigned as either the upper or the lower limit of a new range of possibility, half as large as the previous one. A new trial F is then calculated, in the middle of the new range, and the cycle is repeated. Within ten iterations, the range of possibility is reduced by a factor of a thousand. (210=1024) The true value of F is known to be within this small range. A few more iterations allow F to be found to any desired level of accuracy. This technique is frequently useful.
Greenhouse effect in solar cookers[edit | edit source]
Glass and certain (but not all) "transparent" plastics will allow virtually all the visible and ultraviolet solar radiation striking their surface to pass through, but will absorb most of the infra-red radiation. That is, the visible sunlight comes through the glazing and is directly absorbed by the pot and the air inside, but the thermal energy reradiated by these will not pass directly outward through the glazing. To escape, it must heat the glazing so it radiates the energy outward. This process increases the temperature inside the cooker, as explained below. This greatly aids the cooking. Some facts:
- Sunlight can be considered as having two parts roughly in the same proportion. Broadly the two parts are (1) the visible spectrum light and (2) the infra red rays.
- Both parts can be directly absorbed by the pot and will heat it. However, the colour of the surface of the pot affects how much sunlight is absorbed, and how much is reflected.
- Black paint is good at absorbing visible light, converting its energy into heat. It is also usually good at absorbing infra-red rays. When the pot becomes hot, the black surface re-radiates energy as infra-red rays (black body radiation).
- Glazing materials are transparent to visible light, but fairly opaque to infra-red, especially "far" infra-red, which has wavelengths far from the visible spectrum. When they are hot, these materials radiate infra-red energy, like black bodies.
- Some approximate calculations: Suppose the intensity of sunlight which strikes the glazing is represented by 2x, of which x is visible and the other x is i-r (infra-red). For the cooker to be in equilibrium, the glazing must be hot enough to radiate 2x outward. By symmetry, it also radiates 2x inward, to the pot. The pot therefore receives 2x of i-r, and also x of visible light which has passed unaffected through the glazing. To be in equilibrium, the pot must therefore radiate 3x of i-r outward. (Thus the glazing absorbs x in i-r from the sun, and 3x from the pot. This balances the 2x that it radiates in both directions.) If the glazing were not present, the pot would radiate 2x outward, to be in equilibrium with incoming sunlight. Thus the effect of the glazing is to raise the temperature of the pot so it radiates 50% more than if the glazing were not present. By the fourth-power law of radiation, the pot's temperature measured from absolute zero is raised by about ten percent, since At typical solar-cooking temperatures, this means that the glazing raises the temperature of the pot by about 35 degrees Celsius, or 65 deg. F. This is enough to make a large difference to the effectiveness of the cooker.
- The last paragraph contains many approximations, but its conclusion is roughly right. It is true that glazing greatly improves the effectiveness of a cooker.
- Note that the greenhouse effect depends on the glazing being hot (though not as hot as the cooking pot). If there is a wind which cools the glazing, the pot will be cooled too, even though it is protected from the wind.
- Note also that this calculation makes no mention of reflection. It is often erroneously stated that glazing traps infra-red by reflecting it. This is not true. In fact, gases such as carbon dioxide produce a "greenhouse effect" without reflecting energy at all.
External Links[edit | edit source]
- Green Life Innovators Sunalign page with open source heliostat software, astronomical theory, related programs, etc.
- Wikipedia Solar Cooker page with several calculations about paraboloidal reflectors.
- Wikipedia Rotating Furnace page with proofs that rotating liquid produces paraboloid.
