Get our free book (in Spanish or English) on rainwater now - To Catch the Rain.

# Sun related calculations

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.

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

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.

## Formulae about paraboloidal dishes

### Dimensions of a dish

The dimensions of a symmetrical paraboloidal dish are related by the equation: $\scriptstyle 4FD = R^2$, where $F$ is the focal length, $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: $\scriptstyle P=2F$ and $\scriptstyle Q=\sqrt {P^2+R^2}$, where F and R are defined as above. The linear diameter is then given by: $\scriptstyle S= \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".

The volume of the dish, i.e. the amount of liquid it could hold if the rim were horizontal and the vertex at the bottom, is given by $\scriptstyle V= \frac {\pi R^4} {8 F}$, or equivalently $\scriptstyle V= \frac {1} {2} \pi R^2 D$, where the symbols are defined as above. This 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.

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

### 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 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. (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.