Sun related calculations

From Appropedia
(Redirected from Sun Related Calculations)
Jump to: navigation, search
See also: Sun related texts

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.

Thank you!


Contents

[edit] Code to calculate Equation of Time and Solar Declination

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.

[edit] Computer program

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

[edit] Alternative calculation

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=360/365.24

W is the Earth's mean angular orbital velocity in degrees per day.

A=W*(D+10)

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=A+(360/\pi)*0.0167*sin(W*(D-2))

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: B=A+1.914*sin(W*(D-2)).

C=(A-arctan(tan(B)/cos(23.44)))/180

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=720*(C-nint(C))

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.

[edit] Addendum about solar declination

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:

Declination = - arcsin(sin(23.44)*cos(B))

which is accurate to within a fraction of a degree.

[edit] Program to calculate sun's position and heliostat mirror alignment

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

[edit] Estimating sunrise and sunset data, using analemma

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.

[edit] Earliest and latest sunrise and sunset

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.

[edit] Times of sunrise and sunset

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.

[edit] Azimuths of sunrise and sunset

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.

[edit] Equation of Time, Declination, and Analemma Program

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

[edit] Formulae about paraboloidal dishes

[edit] Dimensions of a dish

The dimensions of a symmetrical paraboloidal dish are related by the equation: \scriptstyle 4FD = R^2, where \scriptstyle F is the focal length, \scriptstyle D is the depth of the dish (measured along the axis of symmetry from the vertex to the plane of the rim), and \scriptstyle R 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: \scriptstyle P=2F (or equivalently \scriptstyle P=\frac{R^2}{2D}) and \scriptstyle Q=\sqrt {P^2+R^2}, where \scriptstyle F, \scriptstyle R, and \scriptstyle D are defined as above. The linear diameter, \scriptstyle L, is then given by: \scriptstyle L= \frac {RQ} {P} + P \ln \left ( \frac {R+Q} {P} \right ), where \scriptstyle \ln(x) means the natural logarithm of  \scriptstyle x , 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 \scriptstyle V= \frac {1} {8 F}\pi R^4 , or equivalently \scriptstyle V=2\pi F D^2, or \scriptstyle V= \frac {1} {2} \pi R^2 D , 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, \scriptstyle \pi R^2 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.

[edit] Paraboloids made by rotating liquids

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 \omega^2 =g, where f is the focal length, \omega 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, \omega in radians per second, and g in metres per second-squared. The angle unit in \omega 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: \scriptstyle fs^2 \approx 447, where f is the focal length in metres, and s is the rotation speed in RPM.

[edit] Focus-balanced reflectors

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".

[edit] Iterative calculations

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.

[edit] Numerical integration

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.

[edit] Binary searches to invert functions

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.

[edit] Greenhouse effect in solar cookers

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 \scriptstyle 1.1^4 \approx 1.5. 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.

[edit] External Links

[edit] See Also