The following text articles were written in connection with computer programs that are listed on the Sun related calculations page of this wiki. These articles and the programs should be read and used together.
The Latitude and Longitude of the Sun
By: David Williams
For solar-energy purposes, it is often important to know, and to be able to predict, where the sun is in the sky at any time. As seen from most latitudes (except in the arctic and antarctic), it rises every day in the eastern part of the sky, and sets in the western part. But, of course, it does some more complicated motions with the passing seasons.
Everyone knows that the sun is further north in the sky in June than it is in December. This variation of the solar latitude, or the sun's "declination" as astronomers like to call it, is responsible for the seasons. Generally, the weather is warmer when the noon-day sun is higher in the sky than when it is lower.
It is less widely known that the sun's longitude, relative to its mean position, also varies with the time of year. The effects of this variation are not as obvious as the effects of the variation of solar declination, but they are significant in some circumstances. For example, there is a difference between "solar time", as shown by a sundial, and "mean time", as shown by a clock. We generally regulate our lives by clocks, so sundials appear to run too fast or too slowly at different times of the year. In early November, a sundial is more than a quarter of an hour "fast", compared with a clock that shows local mean time. By mid-February, it is nearly a quarter of an hour "slow". Not only do sundials lose time during this period, other sun-related events such as sunrise and sunset show the same effect. This is why, for those of us in north-temperate latitudes. the date of the earliest sunset of the winter is several weeks before the date of the latest sunrise. The evenings are already beginning to become brighter when the mornings are still becoming more gloomy.
There are two reasons for this. One is that the earth's orbit around the sun is not exactly circular. Around the time of perihelion, which happens on or about January 3 each year, the earth is closest to the sun and moves slightly faster around its orbit than it does at other times of the year. At aphelion, in early July, the movement is slowest. The speed at which the sun appears to move across the sky each day is essentially the difference between the earth's angular speed of rotation and the angular speed of its orbital motion. When the orbital speed is fastest, at perihelion, the apparent speed of the sun is slowest, so sundials lose time.
The other reason is related to the tilt of the earth's axis. If we imagine that lines of latitude and longitude are marked on the earth's surface, and that its rotation is then stopped, relative to the fixed stars, the path of the sub-solar point, where the sun is vertically overhead, moves around the tilted circle of the ecliptic as the earth moves around its orbit. At the equinoxes, the sub-solar point is at the equator, where the lines of longitude are spaced furthest apart, and is moving at an angle to these lines. At the solstices, the sub-solar point is moving perpendicularly across the lines of longitude, which are closer together than they are at the equator. So the longitude of the sub-solar point changes most rapidly at the solstices and most slowly at the equinoxes. The result when this effect is superimposed on the rotation of the real earth is to make sundials lose time, relative to clocks, near the solstices, and gain time near the equinoxes.
The two reasons act in the same direction during December and January, which are close to both a solstice and perihelion. This is why sundials lose time most rapidly during these months. In June and July, the two effects act in opposite directions, so sundials do not lose time as rapidly then.
The difference between solar time and local mean time is called the Equation of Time. Unfortunately, there is no universally accepted convention as to its sign. The usage I will follow is to regard the Equation of Time as the correction that must be subtracted from the reading of a sundial to get mean time, so it is positive when the sundial is ahead of the clock and vice versa. This is the normal convention. Some writers, however, use the opposite signs. Of course, it does not matter which sign is used, so long as the meaning is clear. In fact, some authors avoid the ambiguity by not using signs at all, but by writing words such as "sundial fast" and "sundial slow" instead.
Understanding the causes of the Equation of Time in a qualitative way is simple. Calculating it quantitatively is a bit more complex.
In this article, I will outline a reasonably accurate way to calculate the sun's longitude and hence the Equation of Time and also the declination (latitude) of the sun on any date in a year. For simplicity, I will ignore all the slow motions of the earth: precessions, nutations, and so on. Over a long period of time, these simplifying approximations will certainly have significant effects. But for this present year, or any year fairly close to it, the results of the calculation should be (and are) quite accurate.
In order to understand the logic behind the calculation, first imagine a planet in a perfectly circular orbit, with its axis perpendicular to the orbital plane, but which is not rotating with respect to the fixed stars. The sub-solar point would move around the equator at a constant speed, taking one year to make a circuit. We can imagine some sort of vehicle that moves so as to remain at the sub-solar point at all times. This vehicle would change its longitude at a uniform speed, so its longitude is directly proportional to time. Let us call this vehicle the "clock car", and imagine that it continues to move around the equator at the same uniform speed even when we make the planet's motions become like those of the earth.
Now let us imagine that the planet's orbit becomes slightly elliptical, but the orbital period remains the same. The angular orbital velocity is no longer quite constant, so the sub-solar point sometimes advances ahead of the clock car, and sometimes falls behind it. The position of the clock car shows "mean time", and the position of the sub-solar point shows "solar time". If the eccentricity of the orbit is small, as is true in the case of the earth which has an orbital eccentricity of only 1.67 percent, the oscillation of the sub-solar point with respect to the clock car is very close to sinusoidal. The amplitude of the variation of the planet-sun distance equals the eccentricity times the mean distance, but the amplitude of the variation of the angular velocity equals twice the eccentricity times the mean velocity. The factor of two comes from the fact that, to conserve angular momentum, the orbital angular velocity is inversely proportional to the square of the sun-planet distance. So we can write equations like:
A = W * T
B = A + 2 * E * SIN(A - SP)
where T shows time, W is the mean orbital angular velocity (so A, which equals W * T, is the longitude of the clock car relative to its longitude when T is zero), E is the orbital eccentricity, B is the longitude of the sub-solar point, and SP is a constant which is zero if we make the value of T zero at perihelion. Otherwise, SP is the angle the planet (the earth) travels around its orbit between the time of perihelion and the time that we choose to use as zero. In this calculation, relating to the earth, I use the December solstice as the zero, which is 12 days before perihelion, and solar days as the units of time, which gives SP the value W * 12.
The longitude of the sub-solar point on our non-rotating planet shows the position of the sun in the planet's sky or, conversely, the position of the planet as seen from the sun. Therefore, the equations above can be used to calculate the angle between the planet-sun vectors at any two times. The difference between the two values of B, at the two times, is this angle.
Let us define a set of Cartesian axes, with the centre of the planet as the origin. The Z axis is perpendicular to the plane of the orbit, the X axis passes through the clock car, and the Y axis is perpendicular to the other two. Let us first calculate the coordinates of the sun in this frame of reference. If we call the distance between the planet and the sun 1 unit, and if the angle between the X axis and the planet-sun vector is B, then the sun's Z is zero, since the Z axis is perpendicular to the orbital plane. The sun's X coordinate is COS(B), and its Y coordinate is SIN(B).
Now let us imagine the X and Z axes and the set of lines of latitude and longitude on the non-rotating planet to be turned by a moderate angle about the Y axis. We have now given our planet an axial tilt. The clock car's position is also rotated, so it continues to move around the equator, which is now tilted with respect to the orbital plane.
This rotation of the frame of reference changes the X and Z coordinates of the sun. (The Y coordinate is not changed, since the rotation is around the Y axis.) The X coordinate was COS(B) and the Z coordinate was zero. After the rotation, through the angle of axial tilt which I will call L (we are already using T for time), the X coordinate becomes COS(B).COS(L) and the Z coordinate becomes COS(B).SIN(L). Actually, since we are talking about the northern-winter solstice, so the angle of tilt is away from the sun, the sun's Z coordinate becomes COS(B).SIN(-L), where L is the positive angle (23.45 degrees) of tilt of the earth's axis.
If we write ST for SIN(-L), we are now in a position to write an expression for the declination of the sun:
SIN(Declination) = COS(B).ST
This is because the sun-planet distance is still 1 unit, and the Z coordinate of the sun is COS(B).ST. Nothing could be simpler! The sun's declination is therefore given by:
Declination = Arcsine(COS(B).ST)
Unfortunately, programming this in a computer language such as QBasic is hindered because the arcsine function is not directly supported. Instead, we have to use a roundabout method using the arctangent function, ATN:
C = COS(B).ST
Declination = ATN(C / SQR(1 - C * C))
We can now also derive an expression for the sun's longitude:
TAN(Longitude) = Y / X = SIN(B) / (COS(B).COS(L)) = TAN(B) / COS(L)
Writing CT for COS(L) (or COS(-L), if we want to use the reverse angle), this becomes:
TAN(Longitude) = TAN(B) / CT
Taking the arctangent of the right hand side of this equation essentially gives an expression for the solar longitude.
We are now very close to an expression for the solar longitude relative to its mean position, which is basically the difference between the longitudes of the sub-solar point and of the clock car. The clock car's longitude is A, so if we calculate the difference between A and the expression for the solar longitude, we get, in BASIC:
Q = A - ATN(TAN(B) / CT)
(I have done the subtraction in the order that gives the conventional sign to the Equation of Time. However, it also makes the positive direction for the longitude westward.)
Now we encounter a little problem. Arctangent is a many-valued function. For any value of the tangent, there are many (an infinite number) of angles, differing from each other by multiples of "half turns", 180 degrees or PI radians. The function ATN produces one of these, but not necessarily the one we want. We may have to add or subtract a multiple of half turns to get the correct angle. Fortunately, in the case of the earth, which has an orbit that is nearly circular, the sun's longitude relative to its mean value is always small. The difference is always less than 0.03 of a half turn. So, if we first divide Q, in the last equation, by PI, to get its value in half turns, and then subtract the nearest integer, we will get an expression for the sun's longitude relative to its mean position with no ambiguity:
C = (A - ATN(TAN(B) / CT)) / PI
Relative_Longitude = C - INT(C + 0.5)
It turns out that the integer that is subtracted is 0, 1, or 2, depending on the time of year.
This expression for relative longitude is in half turns. To get it back in radians, we just have to multiply by PI.
SL = PI * (C - INT(C + 0.5))
Alternatively, we can get the Equation of Time, in minutes, by multiplying by 720 instead of PI, since the earth takes 720 minutes (12 hours) to rotate a half turn relative to the sun.
ET = 720 * (C - INT(C + 0.5))
As the saying goes: The proof of the pudding is in the eating. Theory must always be tested against reality. I wrote a little program that, on 36 days of the year (the 1st, 11th, and 21st days of every month), compares the values of the solar declination and the Equation of Time, calculated as shown above, with values that I got from a published table. (See Reference 1) The agreement between the calculated and published values is striking. The Root Mean Square difference between the calculated and published values of the Equation of Time comes to only 3.7 seconds of time! The greatest difference is 6.0 seconds, on July 11. Compare these with the approximately 130 seconds that the sun takes to move a distance equal to its own diameter in the sky. In the case of the Solar Declination, the R.M.S. difference is 4.7 arc-minutes (60ths of a degree), and the largest difference, on April 11, is 8.8 arc-minutes. At that time of year, the solar declination changes at a rate of about 24 arc-minutes per day. Also, the angular diameter of the sun as seen from the earth is about 32 arc-minutes. The discrepancies between the calculated and published values of the solar declination are only small fractions of either of these.
In my program "SunAlign", which calculates the position of the sun in the sky and the alignment of a heliostat mirror, I used the expressions derived above for the declination and longitude of the sun. These are somewhat better than the approximations that are used in some other programs. For example, some assume that the solar declination varies sinusoidally with time. Others approximate the Equation of Time as the sum of two sine waves, one with a period of a year and the other with a period of half a year. Both of these approximations produce considerable errors.
However, I should not pretend that the expressions derived above are precisely accurate. In addition to the slow variations of the earth's motions mentioned earlier, they do not take account of the non-spherical shape of the earth, the perturbations of its orbit by the moon and other planets, or the effect of parallax when the sun is observed from different points on the earth. Also, the correction for the eccentricity (elliptical shape) of the earth's orbit is only a first-order one. The expressions I have derived are therefore not exact. However, for the purposes of solar energy, they should be amply good enough.
Reference 1
There are many published tables of the Equation of Time and solar declination. The ones I used are from the book:
Sundials, Their Theory and Construction
by Albert E. Waugh
Dover Publications, New York, 1973
ISBN 0-486-22947-5
SunAlign Explained
by David Williams
Computerized Sun Trackers and Heliostats
There are many different designs for machines that track the sun in the sky or move mirrors to reflect sunlight in a constant direction. Some are purely mechanical, using clock mechanisms, others depend on sensors to detect the light or heat of the sun, and yet others use computers to do almost all the work, calculating where the sun is in the sky from astronomical theory, and sending control signals to motors that move the tracker or heliostat mirror. Although the mechanical and sensor-driven ones are conceptually simple, they are often more difficult to make than computerized ones, and have fewer capabilities.
When I built my first heliostat, in about 1987, I decided to make it a computer-controlled one. At that time, BASIC was the language in which all small computers were programmed, so that was the language in which I wrote the software for my machine. The program "SunAlign", whose functioning this article describes, is derived from that heliostat program, and is therefore still in BASIC. Small "microcontroller" computers still use this language, so this program may still be useful in this form.
"SunAlign" is not a complete sun tracker or heliostat program. It has no code for controlling the motors of the machine. But it does do all the astronomical and trigonometrical calculations that are needed.
The program starts with some introductory comments:
_________________________________________________________________
' 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 4Y6
' Canada
' Initially dated 2007 Jul 07
' This version 2008 Nov 03
' All angles in radians except in i/o routines DegIn and DegOut
_________________________________________________________________
The apostrophes at the starts of the lines are used in BASIC to make the computer ignore them during program execution.
Next comes a list of subroutines and functions:
________________________________________________
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)
_________________________________________________
These are tools that the program uses, so I will describe them soon, even though in the real program they are placed at the end of the program. However, there are a few more lines in the main module which we should consider first:
__________________________________________________________________
CONST PY = 3.1415926536# ' "PI" not assignable in some BASICs
CONST PY2 = PY + PY
CONST DR = 180 / PY ' degree / radian factor
_________________________________________________________________
In QBasic, CONSTants are global, so these quantities are accessible within subroutines and functions, as well as in the main module. The strange spelling "PY" is used because some implementations of BASIC require PI (or its equivalent) to be initialized by a value assignment before it is used, but other BASICs have PI as a reserved name, and any attempt to assign a value to it causes an error. Using an unusual spelling allows the program to work in both kinds of BASIC.
The first subroutine to be used is DegIn:
________________________________________________________
SUB DegIn (P$, X)
' Input angle in degrees and convert to radians
PRINT P$;
INPUT N
X = N / DR
END SUB
___________________________________________________________
In QBasic, quantities are passed to and from subroutines "by reference". This means that they can go either way, from the main program to the subroutine or from the subroutine to the main program. DegIn receives a string called P$, which is a prompt to the user. It prints the prompt on the screen, then allows the user to input a number, which is actually the number of degrees in an angle. This number is divided by the constant DR, which is the number of degrees in a radian. The number X is therefore the number of radians in the angle, and this is passed back to the main program.
The subroutine DegOut goes the other way:
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
Again, it receives a prompt string, but this time it also receives a number of radians, X. X is multiplied by DR, to get the angle in degrees, and then some string operations are done to produce a neat output with the angle rounded to the nearest tenth of a degree. Note that any angle greater than 359.95 degrees is converted to 0.0, not to 360.0. North is azimuth zero. The subroutine then prints the output, starting with the prompt.
QBasic does all trigonometrical operations using angles in radians. Therefore, within "SunAlign", all angles are expressed in radians. DegIn and DegOut interface the radians within the program to the degrees with which the user inputs angles and reads them as outputs.
Several operations within the program involve switching three-dimensional co-ordinates between a polar notation, with angles of azimuth and elevation, and a Cartesian one, with X, Y and Z co-ordinates representing the direction of a vector in space. Two subroutines, P2C (Polar to Cartesian) and C2P (Cartesian to Polar), do the conversions.
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
In simple situations, azimuths are the same as true compass bearings, measured clockwise from North. The X-, Y- and Z-axes in the Cartesian notation are directed to the East, North, and upward, respectively. P2C takes an angle of azimuth and an angle of elevation and converts them to equivalent values of X, Y and Z. The resulting vector is of unit length, so the sum of the squares of X, Y and Z is one. X, Y and Z are passed back to the main module.
Note that the "principal axis", the one around which azimuths are measured, is the Z-axis. It is the one that corresponds to the third of the Cartesian co-ordinates in the parameter list that communicates quantities betwen the subroutine and the main module.
Before we consider the sub C2P, which converts in the opposite direction, we should look at the function Ang, which C2P uses.
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
Imagine a point (X,Y) on a plane. The X-axis goes to the right (conventionally) from the origin. The line joining the origin to (X,Y) makes an angle with the X-axis which Ang calculates. SGN, in BASIC, is a number representing the sign of the argument. SGN is 1, 0, or -1, depending on whether the argument is positive, zero, or negative, respectively. ATN is the arc-tangent function, which is always in the range -pi/2 < ATN < pi/2.
Ang is in the range -pi/2 <= Ang < 3pi/2. This is the complete circle.
Some computer languages have built-in functions such as Atan2(X,Y) which are essentially the same as Ang, but which may produce angles in different ranges.
Now we can look at C2P:
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
SQR means square-root.
Ang is used twice, first to calculate the angle of elevation, and then the azimuth. Ang can produce a negative output, down to -pi/2, and this is acceptable as an angle of elevation. The direction in which the mirror must be aimed, for example, can be below the horizontal. But negative azimuths are not used. If Ang produces a negative value for the azimuth, the constant PY2 (2pi) is added to produce a positive equivalent. AZ is therefore in the range 0 <= AZ < 2pi.
Having looked at all the subroutines, we can now return to the main module of the program. After the CONSTants, which we have already seen, it defines some more numbers which are used in the calculations.
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
Perihelion is the date (January 3) when the earth is closest to the sun on its elliptical orbit. The rest of the above definitions should be self-explanatory. Note that since the program uses the northern winter solstice as its zero-point, when the earth's north pole is inclined away from the sun, the value that is used for the axial tilt (or obliquity, in astronomical jargon) is negative.
The next few lines just print explanatory comments.
CLS
PRINT "SunAlign"
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 "Azimuths in degrees (0 to 359.9) clockwise from True North."
PRINT "0 = North, 90 = East, 180 = South, 270 = West."
PRINT "Angles of elevation in degrees upward from horizontal."
PRINT "Negative elevations (downward) allowed for target and mirror."
Then we come to the "Menu", where the user chooses what kind of calculation to do.
Menu:
PRINT "1. Calculate sun's position"
PRINT "2. Calculate mirror orientation"
PRINT "3. Calculate both"
PRINT "4. Quit program"
PRINT "Which? (1 - 4)";
DO
S% = VAL(INKEY$)
LOOP UNTIL S% >= 1 AND S% <= 4
PRINT S%
IF S% = 4 THEN END
The user's choice produces an integer, S%, which must be between 1 and 4, inclusive. However, if it is 4, the program immediately ends, so only the possibilities 1, 2 and 3 remain. S% will be used later in the program to decide what is done.
The next few lines ask the user for various parameters, and lets him input them.
' Note: For brevity, no error checks on user inputs
PRINT "Use negative numbers for directions opposite to those shown."
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%
The Subroutine DegIn, which we have already looked at, is used for the latitude and longitude. LT and LG are therefore in radians, even though the user input them in degrees. The other quantities are simple inputs. Note that the month and day are integers, denoted by the "%" suffix, but all the others can include decimal fractions. There are some places, for example, where the time-zone is not an integer.
Now we come to some astronomical calculations in which the latitude and longitude of the sun are calculated. I won't go into the astronomical theory behind these calculations, but I will summarize what they do.
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
INT rounds a number downward to an integer. MOD does an integer division and discards the quotient, leaving the remainder.
Astronomers call the latitude of the sun, or any other celestial object, its "declination". Note that in the declination calculation a line uses a roundabout method with ATN to calculate the arc-sine of the number C. Arc-sine is not a directly-available function in QBasic. In languages where it is available, the calculation could be done more simply.
The quantity SL is the longitude of the sun relative to the longitude it would have if its motion were exactly uniform, which in reality it is not.
Let us imagine someone who starts from our heliostat or sun tracker and moves due north, along a line of longitude, until he is very close to, but not exactly at, the North Pole. Since he is not quite at the pole, he still has a meaningful longitude, which is the same as that of our machine, and there are still directions that, for him, are North, South, etc.. But he is so close to the pole that the North Celestial Pole is essentially directly overhead for him, and the celestial equator runs around his horizon. So the angle of elevation of the sun in the sky is the same as its declination. The angle LD, which was calculated above, is the difference between the sun's actual longitude (which depends on the time of day) and the direction which our near-polar observer sees as "north". For him, therefore, LD is the compass bearing, or azimuth, of the sun.
Using LD as azimuth and DC as elevation, the next line of the program calls the subroutine P2C to calculate values of X, Y and Z for a unit vector directed toward the sun.
CALL P2C(LD, DC, sX, sY, sZ) ' polar axis (perpend'r to azimuth plane)
The axis around which the azimuth is measured is, of course, directed toward the celestial north pole.
Now suppose that this observer moves south until he is again at the position of our machine. In order for the Z-axis to remain pointing vertically upward, the co-ordinate system must rotate around the X-axis by an angle equal to the co-latitude of the machine, i.e. pi/2 minus its latitude.
The co-latitude, CL, was calculated earlier. The rotation is done in two lines of the program.
CALL C2P(sY, sZ, sX, sAZ, sEL) ' horizontal east-west axis
CALL P2C(sAZ + CL, sEL, sY, sZ, sX) ' rotate by co-latitude
The first of these lines converts the Cartesian co-ordinates back to polar ones, but it is done with the value of X in the third position of the parameter list. This makes the X-axis the principal one in the polar co-ordinates. The azimuth, calculated by that line, is direction of the sun projected onto the Y-Z plane.
The second line converts back to Cartesian notation, but CL is added to the azimuth angle first. This effectively rotates the Y and Z axes around the X one by the angle CL, The Z-axis therefore points vertically upward from the position of our machine.
Since Z is measured upward, a negative value indicates that the sun is below the horizon. This possibility is handled by the next few lines of the program.
IF sZ < 0 THEN
BEEP
PRINT "Sun Below Horizon"
GOTO NewCalc
END IF
The GOTO makes program execution skip over all the remaining calculations, since these are obviously inappropriate if the sun is not visible.
Now we get to the point where the user's initial choices, made in the main menu, become relevant. The value of the integer S% is 1 or 3 if he wants the position of the sun in the sky to be calculated and displayed. S% is 2 or 3 if he wants the alignment of a heliostat mirror to be calculated and displayed.
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
END IF
We already know the X, Y and Z of the vector pointed at the sun. Calculating its angles of azimuth and elevation is just a matter of calling the subroutine C2P. The resulting angles are output to the user as degrees, using the subroutine DegOut, which we looked at earlier.
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
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
END IF
Before the orientation of the mirror can be calculated, the direction of the target, at which light has to be reflected, is input with a couple of calls to DegIn. P2C is then used to find the X, Y and Z of a unit vector directed toward the target.
The direction in which the mirror must be aimed is the bisector of the angle between the directions of the sun and the target. Since we have unit (i.e. equal) vectors directed toward the sun and target, we can calculate the direction of the angle bisector simply by adding the two Xes, the two Ys, and the two Zs. This is done in the CALL C2P line, which then converts to azimuth and elevation. These are then output with DegOut.
At this point, everything important is finished. The remaining few lines of the program just tidy up.
NewCalc:
PRINT "New Calculation"
GOTO Menu
The label NewCalc is the point to which execution jumps if the sun is below the horizon. The final GOTO Menu sends execution back to the main menu, so the user can choose whether to do another calculation or quit the program.
I should point out that SunAlign does make some minor approximations. It takes no account of atmospheric refraction which can change the sun's apparent elevation, but only by a fraction of a degree, and only when the sun is very close to the horizon. Since sunlight is weak when the sun is low, this has little effect for solar energy purposes. Also, the program takes no account of leap years, rounds dates to the nearest day, and a few other small things. These make SunAlign's results slightly approximate, but only by a small fraction of a degree. For solar energy purposes, it is good enough.
I hope people find it useful!
David Williams
Toronto, Canada. 2009 01 23
