Sun related calculations: Difference between revisions

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!

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: [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.

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.

${\displaystyle W=360/365.24}$

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

${\displaystyle 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.

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

${\displaystyle 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⁻¹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:

${\displaystyle 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.

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:

${\displaystyle Declination=-arcsin(sin(23.44)*cos(B))}$

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 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 a specified target. It is accurate to within 0.25 degree.

' 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

Formulae about paraboloidal dishes

Dimensions of a dish

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

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

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.

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 S, 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 S 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 S. Depending on whether the calculated S 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¹⁰=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.

External Links

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

See Also

• Solar photovoltaic software
• Focus-balanced paraboloid

Solar menu

Main topics	Solar energy Solar power Photovoltaics (Research) Concentrated solar power Solar thermal energy Photovoltaic solar thermal Passive solar
Technologies	Solar cooker (Projects) Solar hot water system (Projects) Solar still (Projects) Solar dehydrator (Projects) Greenhouse (Projects) Solar water disinfection (Projects) Photovoltaic systems Photovoltaic devices