Editing Sun related calculations

{{Solar menu}}

{{Sun related calculations}}

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 ==

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 ===

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

Anyone who wants further information should follow the following link: [[http://www.green-life-innovators.org/tiki-index.php?page=Sunalign Link]]
The program that includes the routine is ETIMSDEC. There are instructions how to run it. The article titled "The Latitude and Longitude of the Sun", which I wrote several years ago, describes the astronomical logic behind the routine.

[[User:DOwenWilliams|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 ===

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.

:<math>W=360/365.24</math>

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

:<math>A=W*(D+10)</math>

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.

:<math>B=A+(360/\pi)*0.0167*sin(W*(D-2))</math>

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

:<math>C=(A-arctan(tan(B)/cos(23.44)))/180</math>

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<sup>−1</sup>x'') 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:

:<math>EoT=720*(C-nint(C))</math>

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 ====

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:

:<math>Declination = - arcsin(sin(23.44)*cos(B))</math>

which is accurate to within a fraction of a degree.

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

[[User:DOwenWilliams|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 ==

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 ===

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 ===

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&nbsp;am and 6&nbsp;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&nbsp;am and the time of sunrise on the given date. Measuring the length of this equatorial segment therefore gives the difference between 6&nbsp;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&nbsp;x&nbsp;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&nbsp;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&nbsp;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 ===

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 ==

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 ==

=== Dimensions of a dish ===

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

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

=== Focus-balanced reflectors ===

{{main|Focus-balanced paraboloid}}

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 ==

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 ===

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 ===

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 [http://en.wikipedia.org/wiki/Bisection_method 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. (2<sup>10</sup>=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 ==

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 <math>\scriptstyle 1.1^4 \approx 1.5.</math> At typical solar-cooking temperatures, this means that the glazing raises the temperature of the pot by about 35&nbsp;degrees Celsius, or 65&nbsp;deg.&nbsp;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 ==

* [http://web.archive.org/web/20160331085234/http://green-life-innovators.org/tiki-index.php?page=Sunalign Green Life Innovators Sunalign page] with open source heliostat software, astronomical theory, related programs, etc.
* [http://en.wikipedia.org/wiki/Solar_cooker#Cookers_with_paraboloidal_reflectors Wikipedia Solar Cooker page] with several calculations about paraboloidal reflectors.
* [http://en.wikipedia.org/wiki/Rotating_furnace#Explanation_of_the_equilibrium Wikipedia Rotating Furnace page] with proofs that rotating liquid produces paraboloid.

== See also ==

* [[Solar photovoltaic software]]
* [[Focus-balanced paraboloid]]
* [[Sun related texts]]

{{Page data
| keywords = Solar, Solar energy, solar cooking, parabolic solar cookers, knowledge, mathematics
| sdg = SDG07 Affordable and clean energy
}}

[[Category:Solar energy]]
[[Category:Solar cooking]]
[[Category:Knowledge]]
[[Category:Mathematics]]