NOAA Technical Memorandum ERL ARL-94 

THE NOAA SOLAR EPHEMERIS PROGRAM 

Albion D. Taylor 

Air Resources Laboratories 

Silver Spring, Maryland 

January 1981 
NOAA 'Technical Memorandum ERL ARL-94 

THE NOAA SOLAR EPHEMERlS PROGRAM 

Albion D. Taylor 

Air Resources Laboratories 

Silver Spring, Maryland 

January 1981 
NOTICE

The Environmental Research Laboratories do not approve, 

recommend, or endorse any proprietary product or 

proprietary material mentioned in this publication. No 

reference shall be made to the Environmental Research 

Laboratories or to this publication furnished by the 

Environmental Research Laboratories in any advertising or 

sales promotion which would indicate or imply that the 

Environmental Research Laboratories approve, recommend, or 

endorse any proprietary product or proprietary material 

mentioned herein, or which has as its purpose an intent to 

cause directly or indirectly the advertised product to be 

used or purchased because of this Environmental Research 

Laboratories publication. 
Abstract

A system of FORTRAN language computer programs is presented

which have the ability to locate the sun at arbitrary times. On demand,

the programs will return the distance and direction to the sun, either as

seen by an observer at an arbitrary location on the Earth, or in a standard astronomic coordinate system. For one century before or after the

year 1960, the program is expected to have an accuracy of 30 seconds

of arc (2 seconds of time) in angular position, and 7 105 A.U. in

distance. A non-standard algorithm is used which minimizes the number

of trigonometric evaluations involved in the computations.

1
The NOAA Solar Ephemeris Program

Albion D. Taylor

National Oceanic and Atmospheric Administration

Air Resources Laboratories

Silver Spring, MD

January 1981

Contents

1 Introduction 3

2 Use of the Solar Ephemeris Subroutines 3

3 Astronomical Terminology and Coordinate Systems 5

4 Computation Methods for the NOAA Solar Ephemeris 11

5 References 16

A Program Listings 17

A.1 SOLEFM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

A.2 SOLTIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

A.3 EQ2AZM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

A.4 JULHR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

A.5 DATEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

B The Calendar Routines JULHR and DATEX 22

The Environmental Research Laboratories do not approve, recommend, or endorse any proprietary product or proprietary material

mentioned in this publication. No reference shall be made to the Environmental Research Laboratories or to this publication furnished

by the Environmental Research Laboratories in any advertising or

sales promotion which would indicate or imply that the Environmental Research Laboratories approve, recommend, or endorse any

proprietary product or proprietary material mentioned herein, or

which has as its purpose an intent to cause directly or indirectly the

advertised product to be used or purchased because of this Environmental Research Laboratories publication.

2
1 Introduction

Recently, the Air Resources Laboratories has been working on some projects

connected with the development of a climatology for the ináux of solar energy

at various stations in the U. S. and elsewhere. In connection with these projects,

a computer algorithm was required for locating the sun in position and distance

for arbitrary times and at arbitrary locations on the Earthís surface.

Existing algorithms, such as [Woolf, 1968], do not appear to provide the

distance to the sun. Since the Earth is over 3% closer to the sun in January

than in July, it receives more than 6% more energy ináux then, and this is not

negligible for solar energy considerations.

Many published algorithms ignore the leap year cycle and treat each year as

any other year. This practice simpliÖes the calculation, at a cost in the precision

with which the sun can be located. Because the Earth takes 365.25 days to circle

the sun, the sun has moved an extra quarter day between noon of January 1

of an ordinary year to the same time next year; this produces a discrepancy in

solar time of about half a minute of time, and a discrepancy in declination of

about half a degree. Although slight, these discrepancies were troublesome in

the context of the solar energy climatology.

For these reasons, we have written a system of FORTRAN callable subroutines which locate the solar distance and position for any arbitrary date and

time, keeping proper account of the leap year cycle. We have endeavored to

calculate these positions as closely as practicable. In this context, this is apparently 30 seconds of angle or 2 seconds of time for the angular positions, and

7105 AU for the distance to the sun. Output from our subroutines will di§er

from the published ephemerides by no more than these amounts for at least a

century before and after the year 1966.

Higher accuracy demands consideration of the perturbations caused by the

other planets as well as parallactic changes due to the position on the Earth,

and involve excessive complications. By contrast, the diameter of the sun is

about 9 103 AU and its apparent diameter as seen from Earth is about half

a degree; the point being located is actually the center of the sun.

2 Use of the Solar Ephemeris Subroutines

A listing of the subroutines is found in Appendix A; to use them, the programmer attaches them to his main program and places the statements

JHR = JULHR(KYR, KMON, KDAY, KHOUR)

CALL SOLEFM(JHR, DMIN, RAAPP, DECL, RADVEC, EQTIM)

in the main program.

The INTEGER variables KYR, KMON, KDAY, and KHOUR and the REAL

variable DMIN specify the time, in Greenwich Mean Time (GMT) for which

the Solar position is desired. The year KYR may be either a 4-digit or 2-digit

number (e.g. 1966 or 66), the month KMON may be an integer from 1 to 12,

3
and the day KDAY and the hour KHOUR from 1 to 31 or 0 to 23, respectively.

DMIN speciÖes the minutes and decimal fraction following KHOUR.

Programmers concerned with local standard or daylight time may modify

this sequence for convenience. For example, Eastern Standard Time lags GMT

by 5 hours; a programmer may substitute the call

JHR = JULHR(KYR, KMON, KDAY, KHOUR+5)

where now KYR, KMON, KDAY and KHOUR refer to EST. The function

JULHR will adjust the day, the month, and the year if necessary and compress

this information in the variable JHR for passing on to SOLEFM. For further

information on JULHR, see Appendix B.

If less precision is required, the information on hour, minute, and year may be

replaced by standard values for KYR, KHOUR, and DMIN. Suitable standard

values are 1978, 12, and 0.0; i.e. half through a leap year cycle, at noon GMT.

Values for all other variables will be returned by the subroutine. These

variables and their meaning are as follows:

RAAPP - The Apparent Right Ascension. The angle, in degrees, of the sunís

projection on the celestial equator, as measured from the Vernal Equinox

and corrected for aberration. This is analogous to a pointís longitude on

the Earthís surface.

DECL - The Apparent Declination. The angle, in degrees, of the sun above

(below, if negative) the celestial equator, as corrected for aberration. This

is analogous to a terrestrial pointís latitude.

RADVEC - The radius vector. The distance, in Astronomical Units (AU),

from the Earth to the sun.

EQTIM - The Equation of Time. This is a correction1term representing, in

minutes of time, the di§erence in right ascension between the actual sun

and a Öctitious body called the mean sun. The mean sun moves along the

celestial equator at a uniform rate which is equal, on the average, to that

of the true sun. The equation of time must be added to the Mean Solar

Time (time determined by the mean sun) to obtain the True Solar Time.

These values specify the position of the sun in geocentric, equatorial coordinates. This coordinate system is independent of the position of the observer

on the Earth, and in particular does not specify where the sun is in relation to

the observerís horizon.

Solar energy studies, of course, depend on the relative position of the sun

above the horizon, and a conversion must be made to the site of the observer.

1 The name derives from the centuries old practice among astronomers of using the term

îequationî to mean a correction to be applied to a simple approximation, such as the position

of the mean sun, to get a more accurate expression (the true sun). The term has no particular

relation to the modern mathematical use of the term îequationî.

4
Two subroutines, SOLTIM and EQ2AZM are provided to perform this conversion. They use the information returned by SOLEFM and JULHR and may be

invoked with the statements

CALL SOLTIM( JHR, DMIN, EQTIM, XLONG, STT, STM)

CALL EQ2AZM( DECL, STT, XLAT, ELEV, AZIM)

where XLAT and XLONG are the latitude and longitude of the required site.

Latitude and longitude must be supplied in degrees and decimal fraction, with

North and West positive, South and East negative.

Subroutine SOLTIM returns the following new variables:

STT - True Solar Time. The proportional angle traversed by the sun on its

apparent daily course across the sky, measured in hours and fraction of

time (at 15  per hour, 360 in 24 hours), adjusted so that the crossing of

the meridian (great circle from pole through the zenith) takes place at 12

hours (noon). This is the time (neglecting refraction) which would register

on an ideal sundial.

STM - Mean Solar Time. The value that Solar Time would have if

measured by the mean sun. This will lead or lag civil time by a Öxed

amount.

Subroutine EQ2AZM returns values for the following new variables:

ELEV - Elevation Angle. The angle between the horizon and the center of

the sun, measured in degrees and decimal fraction. Positive values are

returned if above the horizon, negative if below; no correction for refraction

is made.

AZIM - The Azimuth Angle. The angle, measured from True North, eastward

to the horizontal projection of the sun, measured in degrees and fraction.

3 Astronomical Terminology and Coordinate Systems

Before describing the techniques of computation in our subroutines, it is important to acquaint the reader with certain facts and deÖnitions used by astronomers. Only those which have a direct bearing on our program will be

given; for further information, see e.g. [Smart, 1944] or the Explanatory Supplement to the American Ephemeris [USNO, 1977].

If perturbations by other bodies are neglected, the path taken by the Earth

around the sun is an ellipse, and the center of mass of the Earth-sun pair is

located in one focus. This is Keplerís Örst law. To our speciÖed accuracy, the

center of the sun is indistinguishable from that mass-center, and we place the

sun S at one focus in Figure 1, whichdiagrams the orbit of the Earth. For

clarity, the eccentricity of the orbit has been increased by a factor of Öve, and

5
Figure 1: Elliptical path of the Earth about the Sun. Eccentricity exaggerated

for clarity.

the divisions between months adjusted according to Keplerís second law. Even

with this exaggeration, the ellipticity of the orbit is scarcely discernible.

The nearest and furthest points from the sun on the ellipse are called perihelion and aphelion respectively. If the Earth is located at the point E, then

the angle  = PSE is called the anomaly and deÖnes the location of the Earth

on its elliptical path. The ratio e of the distance CS between the center C and

one of the foci S to the distance CP between the center C and either perihelion

or aphelion is called the eccentricity of the ellipse. The eccentricity speciÖes the

exact shape of the ellipse.

The mean distance from sun to Earth is deÖned to be the distance PC = AC

= half the sum of PS+AS, and this serves to deÖne the length of an astronomic

unit (AU). For historic reasons, 1 AU is actually slightly less (by 30 x 109)

than the mean distance, but this is not signiÖcant for present purposes.

6
Figure 2: Comparison of the Celestial (ecliptic) and Equatorial Coordinate Systems

An ellipse is a plane curve, and so from the viewpoint of an observer on the

Earth, the sun appears to remain in a plane passing through the Earth and

called the plane of the ecliptic. The entire sky can be projected onto a celestial

sphere having the observer as center; the intersection of the plane of the ecliptic

with that sphere (see Figure 2) is a great circle called the ecliptic.

Because of the Earthís diurnal rotation, the stars all appear to describe daily

circles centered on two poles on the celestial sphere, which are projections of

the Earthís North and South Poles. Similarly, the Earthís equator is projected

onto a celestial equator.

Because of the orientation of the Earthís axis, the equator and ecliptic meet

at only two points, called the equinoxes. The angle " between the equator and

ecliptic at the equinoxes is called the obliquity of the ecliptic. The apparent

passage of the sun through one of the equinoxes, (the Vernal Equinox V) to

the Northern Hemisphere marks the beginning of spring; through the other (the

Autumnal Equinox F) marks the beginning of Fall.

The ecliptic and the equator each form the basis of a set of spherical coordinates for specifying the positions of points in the sky such as stars or planets.

Given any point Q (see Figure 2), drop an arc of a great circle meeting the

7
Figure 3: Comparison of the Topocentric Equatorial and Azimuthal Coordinate

Systems

ecliptic at a point Qí. The angular size in degrees of the arc QQíis called the

celestial latitude of Q. The celestial longitude is the arc VQí, measured from

the Vernal Equinox to Qíin the direction of motion of the sun, also in degrees.

These form the ecliptic or celestial coordinates. By convention, the latitude of

a celestial object is denoted by  and the longitude by .

Another system of coordinates is the system of equatorial coordinates. An

arc QQîis dropped to the Equator; its angular size is the declination. The arc

VQî from the Vernal equinox to Qî in the direction of motion of the sun is

called the right ascension. By convention, the declination of a celestial object

is denoted by  and its right ascension by .

The process of conversion from one system to the other is simpliÖed by the

use of the same point (the Vernal Equinox) as a zero for both.

The celestial and equatorial coordinate systems are used to locate the position of one celestial point relative to other celestial points. In order to relate

that celestial point to the reference system of an Earthbound observer, we need a

means of specifying coordinates relative to that observer (topocentric systems).

To an observer located at a prescribed latitude and longitude on the Earth,

the celestial sphere appears tilted so the North celestial pole is at an angle

above the horizon equal to his latitude (below, for South latitudes). The point

Z immediately overhead is called the Zenith; the half great circle from the North

Pole through the Zenith and the horizon to the South pole is called the Meridian

8
(see Figure 3). These reference points and circles depend on the observerís

location and are di§erent for observers located elsewhere.

For a given point Q, the angle from the Meridian to Q subtended at the North

pole and increasing Westward is called the Hour Angle. The hour angle and

the declination (conventionally denoted by h and , respectively) form another

system of equatorial spherical coordinates, one which depends on the location

of a speciÖc observer. Indeed, the hour angle of a point as measured by this

observer is less than the hour angle measured simultaneously by an observer

on the Greenwich Meridian in an amount equal to his (West) longitude. The

combination of Hour angle and Declination form a set of topocentric equatorial

coordinates.

The Önal set of spherical coordinates of interest to us is the system of azimuthal coordinates based on the horizon. In this system, the arc dropped from

Q normal to the horizon is in angular measure the altitude a (also called elevation), and the angular measure of the arc along the horizon from true North

eastward is the azimuth A.

Conversions between the celestial ecliptic and equatorial system, and between the topocentric equatorial and the azimuthal coordinates are straightforward matters of spherical trigonometry. The formulae involved are independent of the time and of the observerís position. To go from the celestial to the

topocentric, the formulae do depend on time and position. The simplest case

is between the two equatorial systems, since the declination of a point does not

change in the transition, and the only change is between right ascension and

hour angle.

In making the change between right ascension and hour angle, one needs

the hour angle of the vernal equinox, to which one can add the right ascension

of the point in question. Like the hour angles of all other celestial points, this

changes by 360 in somewhat less than 24 hours, and may be considered as a

celestial îclockî. The value of the hour angle of the Vernal Equinox is called

the sidereal time provided it is measured in time units. It will change with the

longitude, and an observer must subtract his (west) longitude from sidereal time

at Greenwich to Önd sidereal time at his own site. Thus, if an observer knows

his longitude, the sidereal time at Greenwich, and the right ascension of a point,

he may Önd sidereal time at his site by subtracting the longitude, and the hour

angle of the point by adding the right ascension.

Alternatively, the hour angle of the sun is called (true) solar time; if this

is known at Greenwich together with the sunís right ascension, sidereal time

may be obtained by subtracting the right ascension from the solar time and one

may proceed as before. By convention, an imaginary point is considered which

travels at a uniform rate along the equator at the same average speed as the sun

does along the ecliptic; it may be considered as what the sun would do if the

Solar system were simpler. This point is called the mean sun and its hour angle

mean solar time. At Greenwich, mean solar time is approximately universal

time, on which civil time is based, and (to date) di§ers by less than one minute

from Ephemeris time on which the American Ephemerides of the sun and moon

are based; to the accuracy with which we are concerned, all three time scales

9
may be considered as Greenwich Mean Time.

Thus, approximating Ephemeris time with Greenwich Mean Time, the right

ascension of the mean sun is calculated on the basis of the formula deÖning it;

this leads to the sidereal time by subtraction and then to the hour angle of a

point such as the sun when its right ascension is known.

Because of the use of the hour angles in deÖning time scales, they (together

with the right ascension) are usually measured in hours, minutes and seconds

of time, rather than degrees, minutes and seconds of angle. To avoid confusion

between the two types of minute and second, minutes and seconds of time are

usually written with superscript m and s, rather than the single and double

marks. The size of the time-measure quantities is 15 times that of the degree

measure quantities. Thus, 1

h = 15 

,1

m = 150 and 1s = 15".

Ideally, the major features described above such as the equator, ecliptic,

equinoxes, eccentricity and obliquity would be unchanging constants. In actual

fact, they change slowly with time.

Due to the asphericity of the Earth together with the actions of the other

planets, the axis of the Earthís rotation itself traces a small circle through the

sky in about 23,000 years; this process is called precession. Due to precession,

the equator slowly shifts and so do the equinoxes. Since the Vernal equinox is

the origin for the celestial coordinates, the longitude and right ascension of the

îÖxed starsî also shift.

One result of this is that the right ascension of the perihelion (alternately,

the anomaly of the equinox) decreases, and the Earth takes more time to move

from one perihelion to the next than it takes to move from one equinox and

back again. The time it takes to move from perihelion to perihelion (i.e. for the

anomaly to change by 360 ) is called the anomalistic year. The time it takes to

move from equinox to equinox (for the longitude to change by 360 ) is called a

tropical year and is the basis for the ordinary civil calendar since the equinoxes

mark the seasons. Another often used year is the sidereal year, the time it takes

to move 360 relative to the stars.

The perturbations induced on the motion of the Earth by the moon and the

planets are small, but a portion does accumulate. The resulting accumulated,

long term changes in, inter alia, the eccentricity and obliquity of the Earthís

orbit are known as secular variations.

Many parameters, such as anomaly or longitude may be regarded as made

up of two parts; a mean part due to long term motions (which are either uniform

or include secular variations), and another part due to short term áuctuations.

Thus, we may refer to mean anomaly, mean longitude, etc. The total is called

the true value, the di§erence between true and mean is called an equation. For

example, the di§erence between true anomaly and mean anomaly is called the

equation of the center, while the di§erence between true solar time and mean

solar time is the equation of time.

Astronomers specify time in terms of Julian Date, which is the number

of days and decimal fraction since the date January 1, 4713 B.C. (old Julian

calendar). This practice simpliÖes calculating motion between two given dates,

avoiding the irregularity in the lengths of the months and years. The Julian day

10
begins at noon in civil or universal time, for the convenience of astronomers.

4 Computation Methods for the NOAA Solar

Ephemeris

In this section, we give a summary of the theoretical basis on which the SOLEFM

subroutine was written. This section need not be read to understand the use of

the routine, but it can be useful for comparison with other similar routines.

The physical basis for the computations is that of the dynamics of orbits

as established from the works of Kepler through Newcomb. Similar techniques

have been used for centuries to compute the ephemerides (tabulations of forecast

positions of the Sun, moon, and planets). However, whereas the published

ephemerides include the e§ects of perturbations on the Earthís orbit by the

major planets and the moon, we will ignore these terms and consider a simpliÖed,

two body problem.

We shall use the Earth orbital data from Newcombís tables [Newcomb, 1898]

and the Sunrise and Sunset tables [USNO, 1945] and methods of calculation for

two-body problems adapted from [Smart, 1944], [Pollard, 1976] and [Wintner, 1941].

See these references for further details or for information on techniques for correcting for perturbations by other bodies.

From the Sunrise and Sunset tables, and from Newcombís tables, we Önd the

orbital parameters for Earth, together with their daily changes, for the epoch

1966 January 0.75 UT (Julian Day 243 9126.25), i.e. the positions and velocities

of the Earth relative to the sun at (approximately) 1800 GMT December 31,

1965.

The orbital parameters and their values are given in Table 1.

TABLE 1 - Earth Orbital Parameters

1966 JAN 0.75 UT

PARAMETER VALUE DAILY CHANGE

geometric mean longitude 279 .95656 0



.985 647 3463

Earthís mean anomaly 357 .60087 0



.985 600 2614

eccentricity of Earthís orbit 0.016723401 -1.115 x 109

obliquity of ecliptic 23 .44371 -0 .000 000 35626

right ascension of mean sun 279 .95232 0



.985 647 3494

The secular variations in obliquity and eccentricity were computed from

Newcombís tables; all other terms were taken from the Sunrise and Sunset

tables.

According to Newcombís tables, the largest perturbations in angular positions of the sun are caused, in order, by Jupiter, the Moon, Venus, Mars, Saturn,

and Mercury; together they cause variations of up to about 30 seconds of arc

or two seconds of time. For perturbations of the radius vector, the important

bodies are, in order, the Moon, Jupiter, Venus, Mars and Saturn; together they

cause variations of up to 7 x 105 AU. These values, then, constitute the limits

of accuracy of our subroutines as stated earlier. Over a century, the secular

11
variations in obliquity and eccentricity accumulate to more than these limits,

and we have therefore included them. The daily changes also undergo secular

variations, but they are not signiÖcant over a century or less.

To obtain the current mean value of any of the terms listed in Table 1, simply

calculate d, the number of days and decimal fraction elapsed since the stated

epoch, multiply by the daily change, and add to the tabular value. Thus, for

example, the mean anomaly M is given by

M = 357 :60087 + d 0



:9856002614

Routines such as JULHR are invaluable for the calculation of d.

In order to calculate the true anomaly and Solar distance from the mean

anomaly, we recall Keplerís Örst two laws for orbits in a two body system:

1. The path of one body about the other is that of an ellipse and the center

of attraction is in one of the focal points of the ellipse.

2. The radius vector between the two bodies sweeps out areas of the ellipse in

proportion to the time elapsed. (This is the law of conservation of angular

momentum.)

In Figure 4 is diagrammed this elliptical path of the Earth. Suppose the

semi-major axis CP is unity (in fact, very close to one AU), then the distance

SC from one focus (sun) to the center C is e units, e being the eccentricity. The

semi-minor axis is then equal to f =

p

(1  e

2) .

Consider a circle of center C and unit radius circumscribed about the ellipse.

The ellipse may be regarded as a uniform compression of that circle in a scale

factor f in the direction normal to the major axis. Thus, if the Earth is located

at E, and EíEEîis drawn normal to the major axis meeting the axis in Eîand

the circle in Eí, then the ratio of EîE to EíEîis f.

The angle PCEíis known as the eccentric anomaly  and provides a useful

intermediary between the mean anomaly M and the desired true anomaly .

The area of the sector PSEícan be shown to be



2



e sin()

2

(where  is in radian measure). Since the area of the sector PSEíis equal to f

times that of the sector PSE, and thus by Keplerís second law proportional to

the mean .cc 5 anomaly M, we have (in radian measure)

M =   e sin() (1)

This is known as Keplerís equation.

Further, from Figure 4, it can be shown that

r sin() = f sin() (2)

r cos() = cos()  e

12
Figure 4: Relations between the various Anomalies of the Earthís orbit. Eccentricity greatly exaggerated

where r is the radius vector (Earth-Sun distance). From 2, we may also calculate

r = 1  e cos() (3)

Equation 1 may be solved by any of a number of iterative schemes. The one

chosen for these subroutines is Newtonís method. In this case, we set initially

0 = M

and, for k>0,

k = k1 +

M  k1 + e sin(k1)

1  e cos(k1)

until kconverges to the required value .

13
It can be shown that convergence is guaranteed for e < 1:0 (i.e., when the

orbit is an ellipse). For the low eccentricities of the Earth sun pair, the convergence is extremely fast. Indeed, in this case, the Örst iterate 1

is always within

0.5îof the exact solution, well within the accuracy limits we have speciÖed, and

to within this accuracy we may write

 = M +

e sin(M)

1  e cos(M)

(4)

Using the eccentric anomaly , we may Önd the Solar distance r from 3 and

the true anomaly  from 2. The next step is to Önd the location in the various

coordinate systems. The sunís geometric longitude g is found by subtracting

from  the anomaly of the equinox, i.e. the di§erence between the mean anomaly

M and the sunís mean longitude m:

g =   M + m

From this, we subtract a correction for aberration to obtain the apparent longitude.

Aberration is an apparent shift in the position of a celestial object due to the

Önite speed of light and the motion of either the object or the observer. From

the moving Earth, the sun appears shifted slightly in the direction opposite to

the Earthís motion. The amount is small, namely b = 20":50 divided by the

solar distance, but is a large enough fraction of our accuracy limit to be worth

including. Thus, the apparent longitude  is given by

 = g 

b

r

The latitude of the sun, both geometric and apparent, is always zero, since

the Earth-sun line and the terrestrial motion are always (neglecting perturbations) in the plane of the ecliptic. Accordingly, following the formulae for conversion between the ecliptic and equatorial systems given in [USNO, 1977]pp24§,

we have:

sin() = sin(") sin()

cos() = cos()

cos()

(5)

sin() = cos(") sin()

cos()

where  and  denote the declination and right ascension (apparent) of the sun,

respectively, and " denotes the obliquity of the ecliptic.

From the right ascension  we subtract the right ascension of the mean sun

m calculated from Table 1, to obtain the equation of time. This is the last of

the terms returned by Subroutine SOLEFM.

The above algorithm di§ers from the traditional practice of expanding the

expression   M, known as the equation of the center, in terms of a Fourier

14
sine series in the mean anomaly M. Actually, the di§erence between any pair

of anomalies is an odd periodic function of any anomaly and can be expressed

as a sine series of that anomaly. [Smart, 1944] demonstrates how to utilize an

iterative solution of 1 with an assumed series to evaluate the coe¢ cients of  -

M, then to convert that into a series expansion for the equation of the center.

The results are coe¢ cients that are functions of the eccentricity.

The cost of the traditional method is a signiÖcant amount of analytic e§ort

in evaluating the coe¢ cients. There are two major beneÖts. The Örst is the

fact that the coe¢ cients, when calculated for other planets as well, provide a

good starting point for calculating the e§ects of the perturbations induced on

the Earth by those planets. The second is the fact that, all the analytic e§ort

being already expended, such an expression is well suited for hand calculation,

consisting of a series of table lookups, multiplications and addition.

The Örst beneÖt is inapplicable in the present case, since we have elected

to neglect perturbations. As to the second, for computer based algorithms

the computation of trigonometric functions is relatively costly and should be

minimized. The algorithm presented in this paper involves fewer sine evaluations

than any Fourier sine series of comparable accuracy, and there is no need to

reevaluate coe¢ cients for di§ering eccentricities.

Persons interested in accounting for perturbations, or producing ephemerides

for the various planets, are referred to [VanFlandern & Pulkkinen, 1979], where

conventional trigonometric terms are supplied.

The complete account of the solar radiation incident on the ground requires

consideration of atmospheric and meteorological factors as well. These factors

include refraction, attenuation, and climatological e§ects. The present paper

concentrates on the astronomical aspect of the question, and we will not deal

with these other matters.

Refraction is signiÖcant primarily at low solar elevation angles, where the

attenuation is greatest and there is little solar energy available in any case. At

these low angles, it is of the order of 34íof arc and will a§ect the times of sunset

and sunrise by about 2 minutes of time (more at higher latitudes). Refraction

depends on the elevation angle of the sun, as well as temperature gradients and

the altitude of the observer. An analysis may be found in [Smart, 1944].

There are a number of models of the attenuation of solar radiation as it

passes through the atmosphere. In general, they depend on the total air mass

through which the radiation passes, and except for very low elevation angles,

this is approximately inversely proportional to the sine of the elevation angle.

Some models may be found in, e.g.[Hoyt, 1979] and [Reivfeim, 1978].

Climatological e§ects include atmospheric turbidity, cloud cover, and atmospheric moisture. The study of these processes and their a§ects on solar

radiation is the subject of much current research. Some results may be found,

e.g. in [Cotton, 1979] .

15
5 References

References

[Cotton, 1979] Cotton, Gerald F., îARL Models of Global

Solar Radiationî, Appendix VI in SOLMET, Volume 2, Final Report, TD-9724

Hourly Solar Radiation - Surface Meteorological Observations. Environmental Data

and Information Services, NOAA, Asheville,

NC 28801 (1979)

[Hoyt, 1979] Hoyt, Douglas V., îTheoretical Calculations

of the True Solar Noon Atmospheric Transmissionî, Appendix V in SOLMET, Volume

2, Final Report, TD-9724 Hourly Solar Radiation - Surface Meteorological Observations.

Environmental Data and Information Services, NOAA, Asheville, NC 28801 (1979)

[Newcomb, 1898] Newcomb, Simon, îTables of the Motion

of the Earth on its Axis and Around the

Sunî, Astronomical Papers of the American

Ephemeris, Vol VI, part I U.S. Naval Observatory, (1898)

[Pollard, 1976] Pollard, Harry, îCelestial Mechanicsî, Carus

Mathematical Monographs, No. 18, Mathematical Association of America, (1976)

[Reivfeim, 1978] Revfeim, K.J.A. îA Simple Procedure for Estimating Global Daily Radiation on Any Surfaceî, Jour. Appl. Meteor. vol 17, no 8, pp.

1126-1131 (1978)

[Smart, 1944] Smart, W.M., îTextbook on Spherical Astronomyî, Cambridge, (1944)

[USNO, 1945] USNO îTables of Sunrise, Sunset, and

Twilightî, supplement to the American

Ephemeris, 1946, U.S. Naval Observatory,

(1945)

[USNO, 1977] USNO îExplanatory Supplement to the

American Ephemeris 1978î, U.S. Naval Observatory, (1977)

16
[VanFlandern & Pulkkinen, 1979] Van Flandern, T.C. and K.F. Pulkkinen,

îLow-Precision Formulae for Planetary Positionsî, Astrophysical Journal Supplement

Series, vol 41; pp.391-411, (1979)

[Wintner, 1941] Wintner, Aurel, îThe Analytical Foundations of Celestial Mechanicsî, Princeton University Press, (1941)

[Woolf, 1968] Woolf, Harold, îOn the Computation of Solar Elevation Angles and the Determination

of Sunrise and Sunset TimesîNASA TM X1646, NTIS, SpringÖeld, VA (1968)

A Program Listings

A.1 SOLEFM

SOLEFM Locates the Sun in Celestial Coordinates

SUBROUTINE SOLEFM(JHR,DMIN,RAAPP,DECL,RADVEC,EQTIM)

C INPUT TO SUBROUTINE: JHR (JULIAN HOUR FROM SUBROUTINE JULHR),

C DMIN (MINUTES PAST HOUR); ASSUMED VALUES IN EPHEMERIS TIME

C (APPROXIMATELY GREENWICH MEAN TIME).

C RETURNED VALUES: RIGHT ASCENSION AND DECLINATION OF APPARENT

C SUN (AS CORRECTED FOR ABERRATION), DISTANCE TO GEOMETRIC SUN

C (NOT CORRECTED FOR ABERRATION), AND THE EQUATION OF TIME (THE

C HOUR ANGLE OF THE APPARENT SUN MINUS THAT OF THE MEAN SUN, A

C FICTITIOUS POINT MOVING ALONG THE CELESTIAL EQUATOR AT UNIFORM

C RATE AND DEFINING EPHEMERIS TIME).

C

C UNITS: DECLINATION AND RIGHT ASCENSION IN DEGREES AND FRACTION,

C RADIUS VECTOR IN ASTRONOMIC UNITS (AU) AND FRACTION, THE

C EQUATION OF TIME IN MINUTES OF TIME AND FRACTION. (ONE MINUTE

C OF TIME = 15 MINUTES OF ARC.)

DOUBLE PRECISION GML,GML0,GMLPR,ASL,EMA,EMA0,EMAPR

DOUBLE PRECISION EQCENT,SOLDST,OBLQ,OBLQ0,OBLQPR

DOUBLE PRECISION TAU,TAU0,TAUPR,ECCEN,ECCEN0,ECENPR

DOUBLE PRECISION DAYS,RADPDG,ABRCON,ECA,COSECA,SINECA,TANOM

DOUBLE PRECISION F,X0,Y0,X1,Y1,DARSIN,DARG

LOGICAL READY

DATA GML0/279.95656D0/,GMLPR/.985 647 3463D0/

DATA EMA0/357.60087D0/,EMAPR/.985 600 2614D0/

DATA OBLQ0/23.44371D0/,OBLQPR/-3.562 6283D-7/

DATA TAU0/279.952 23D0/,TAUPR/.985 647 3494D0/

17
DATA ECCEN0/1.672 3401D-02/,ECENPR/1.1149D-9/

DATA RADPDG/1.745 329 252D-2/

DATA ABRCON/5.693 333 333D-3/

DATA READY/.FALSE./

DARSIN(DARG)=DATAN2(DARG,

A DSQRT(DMAX1((1D0+DARG)*(1D0-DARG),0D0)))

IF(READY) GO TO 5

JHEPOK=JULHR(1966,1,0,18)

READY=.TRUE.

5 DAYS=(DBLE(FLOAT(JHR-JHEPOK))+DMIN/60D0)/24D0

C DAYS ELAPSED SINCE (OR PRECEDING) EPOCH OF ORBITAL ELEMENTS.

GML=(GML0+DAYS*GMLPR)*RADPDG

C MEAN GEOMETRIC LONGITUDE OF SUN, MEASURED ON ECLIPTIC FROM

C VERNAL EQUINOX.

EMA=(EMA0+DAYS*EMAPR)*RADPDG

C EARTHíS MEAN ANOMALY, MEASURED AT SOLAR CENTER FROM PERIHELION

C TO FIDUCIAL POINT REPRESENTING EARTH, BUT MOVING AT UNIFORM

C SPEED AND TOUCHING EARTH AT APHELION AND PERIHELION.

OBLQ=(OBLQ0+DAYS*OBLQPR)*RADPDG

C OBLIQUITY OF THE ECLIPTIC

TAU=(TAU0+DAYS*TAUPR)*RADPDG

C POSITION OF MEAN SUN ON THE CELESTIAL EQUATOR

ECCEN=ECCEN0+DAYS*ECENPR

F=DSQRT(1D0-ECCEN**2)

C ECCENTRICITY OF THE TERRESTRIAL ORBIT

C ALL THE ABOVE ELEMENTS OF MEAN TERRESTRIAL ORBIT EVALUATED FOR

C 1966 JAN 0.75 UT. (JULIAN DAY 2439126.25). SEE

C SUNRISE AND SUNSET TABLES SUPPLEMENT TO THE AMERICAN

C EPHEMERIS, 1946.

ECA=EMA

COSECA=DCOS(ECA)

SINECA=DSIN(ECA)

ECA=EMA+ECCEN*SINECA/(1D0-ECCEN*COSECA)

C ECA = EARTHíS ECCENTRIC ANOMALY.

COSECA=DCOS(ECA)

SINECA=DSIN(ECA)

TANOM=DATAN2(F*SINECA,COSECA-ECCEN)

C TANOM = EARTHíS TRUE ANOMALY.

SOLDST=1D0-ECCEN*COSECA

RADVEC=SOLDST

EQCENT=TANOM-EMA

C EQCENT = EQUATION OF THE CENTER = TRUE EARTH ANOMALY MINUS

C MEAN EARTH ANOMALY.

ASL=GML+EQCENT-ABRCON*RADPDG/SOLDST

C APPARENT SOLAR LONGITUDE = GEOMETRIC MEAN LONGITUDE PLUS

C EQUATION OF THE CENTER MINUS ABERRATION.

18
Y0=DSIN(ASL)

DECL=DARSIN(DSIN(OBLQ)*Y0)/RADPDG

C DECL = DECLINATION OF THE APPARENT SUN IN DEGREES AND FRACTION;

C RETURNED.

Y0=Y0*DCOS(OBLQ)

X0=DCOS(ASL)

RAAPP=DATAN2(-Y0,-X0)/RADPDG+180D0

C RAAPP = RIGHT ASCENSION OF THE APPARENT SUN IN DEGREES AND

C FRACTION; RETURNED.

Y1=DSIN(TAU)

X1=DCOS(TAU)

EQTIM=4D0*DATAN2(Y1*X0-X1*Y0,Y0*Y1+X0*X1)/RADPDG

C EQUATION OF TIME IN MINUTES OF TIME AND FRACTION; RETURNED.

RETURN

END

A.2 SOLTIM

SOLTIM - Returns True and Mean Solar Time

SUBROUTINE SOLTIM(JHR,DMIN,EQTIM,XLONG,STT,STM)

GMTIM=FLOAT(MOD(JHR,24))+DMIN/60.

STM=GMTIM-XLONG/15.

STM=AMOD(AMOD(STM,24.)+24.,24.)

STT=STM+EQTIM/60.

STT=AMOD(STT+24.,24.)

RETURN

END

A.3 EQ2AZM

EQ2AZM - Converts from Equatorial to Azimuthal Coordinates

SUBROUTINE EQ2AZM(DECL,STT,XLAT,ELEV,AZIM)

DATA RADPDG/1.745329252E-2/,RADPHR/2.617993878E-1/

DECLR=RADPDG*DECL

XLATR=RADPDG*XLAT

STTR=RADPHR*STT

SLAT=SIN(XLATR)

CLAT=COS(XLATR)

19
SDECL=SIN(DECLR)

CDECL=COS(DECLR)

CSTT=COS(STTR)

SSTT=SIN(STTR)

Z=SLAT*SDECL-CDECL*CLAT*CSTT

X=-CDECL*SSTT

Y=-CLAT*SDECL-CDECL*SLAT*CSTT

R=SQRT(X**2+Y**2)

IF(R.EQ.0.) GO TO 20

AZIM=ATAN2(X,Y)/RADPDG+180.

ELEV=ATAN2(Z,R)/RADPDG

RETURN

20 ELEV=SIGN(90.,Z)AZIM=0.

RETURN

END

A.4 JULHR

JULHR - Returns Hours Since Beginning of Julian Calendar

FUNCTION JULHR(MYR,KMO,KDA,KHR)

DIMENSION MONTH(12)

DATA MONTH /0,31,60,91,121,152,182,213,244,274,305,335/

KYR=MYR

IF(MIN0(KYR,99-KYR).GE.0)

KYR=KYR+1900

LMO=KMO

LHR=LMO/12

LMO=LMO-12*LHR

IF(LMO.GT.0) GO TO 10

LMO=LMO+12

LHR=LHR-1

10 KYR=KYR + LHRLHR=KHR+24*(KDA+MONTH(LMO)+366*(KYR-2000))

IF(LHR.LT.1464) GO TO 22

INCR1=(LHR-1464)/8784

GO TO 25

22 INCR1=(LHR-1463)/8784

C INCR1= NUMBER OF YEARS PASSED SINCE 2000 MARCH 1 00Z

C IF NEGATIVE, UNTIL 2000 FEB 29 23Z

25 INCR2=INCR1/4

C INCR2= NUMBER OF YEARS THAT ARE LEAP YEARS

INCR3=INCR2/25

C INCR3= NUMBER OF CENTURIES PASSED

20
INCR4=INCR3/4

C INCR4 = NUMBER OF 4- CENTURY PERIODS

JULHR=LHR+(245 1544 - INCR1 + INCR2 - INCR3 + INCR4) * 24

RETURN

END

A.5 DATEX

DATEX - Converts JULHR Value to Calendar Date

SUBROUTINE DATEX(JULHR,KYR,KMO,KDA,KHR,KNAMMO,KDAYWK)

INTEGER WEEK

DIMENSION MONTH(12),NAMEMO(12),WEEK(7)

DATA MONTH /0,31,60,91,121,152,182,213,244,274,305,335/

DATA NAMEMO/3HJAN,3HFEB,3HMAR,3HAPR,3HMAY,3HJUN,3HJUL,

1 3HAUG,3HSEP,3HOCT,3HNOV,3HDEC/

DATA WEEK/3HSUN,3HMON,3HTUE,3HWED,3HTHU,3HFRI,3HSAT/

MDA=JULHR/24

KHR=JULHR-MDA*24

KDAWK=MOD(MDA+1,7)

KDAYWK=WEEK(KDAWK+1)

MDA=MDA-245 1605

C NUMBER OF DAYS SINCE (NEGATIVE, UNTIL) 2000 MAR 1

IF(MDA.LT.0) GO TO 5

N400= MDA/146 097

GO TO 10

5 N400=(MDA+1)/146 097 -1

C NUMBER OF 400-YEAR PERIODS SINCE DAY 2000 MAR 1 (IF NEGATIVE,

C UNTIL 2000 MAR 1)

10 MDA=MDA-146 097*N400

N100=MDA/36 524

C N100= NUMBER OF CENTURIES SINCE LAST 400-YR PERIOD

MDA=MDA-36 524 * N100

NO4=MDA /1461

MDA=MDA- 1461 *NO4

NYR= MDA/365

MDA=MDA-365*NYR+61

KYR=2000+NYR+4*(NO4+25*(N100+4*N400))

IF(MAX0(NYR,N100).LT.4) GO TO 20

C LEAP DAY

KMO=2

KDA=29

KNAMMO=NAMEMO(2)

21
RETURN

20 IF(MDA.LE.366) GO TO 30

KYR=KYR+1

MDA=MDA-366

30 DO 40 K=2,12

IF(MONTH(K).GE.MDA) GO TO 50

40 CONTINUE

K = 13

50 KMO = K-1

KDA=MDA-MONTH(KMO)

KNAMMO=NAMEMO(KMO)

RETURN

END

B The Calendar Routines JULHR and DATEX

In the automatic processing of archived data covering a period of years, a common source of problems to programmers lies in the nature of our calendar,

particularly the irregular number of days per month. This presents especial

di¢ culties if an analysis calls for data points uniformly spaced in time. Again,

if one wishes to test a sequence of points for missing data, it makes a di§erence

whether the absence of an entry for the 31st is due to the lack of data or the

lack of a day. Even when a programmer has carefully tested each month for its

expected number of days, he may come to Önd on Feb. 29 that he forgot about

leap year.

Historically, this irregularity has been of annoyance to astronomers, who

may well need to know how many days it is between, say, Jan. 1, 1900 and

Dec. 21, 1960. Although the arithmetic involved is not sophisticated, it can be

quite tedious, and it is not uncommon to make errors of perhaps several days.

If many such calculations are required, the calculation easily goes from tedious

to costly.

The solution to this problem was provided in 1582 by the Italian Protestant

scholar, Joseph Scaliger. Basically, his system was to assign to each day a

number, the Julian Date, which was the number of days elapsed since Jan. 1,

4713 B.C., this date being arbitrarily chosen as being early enough to precede all

historical events and precisely observed and recorded astronomical phenomena.

The name Julian was applied in honor of his father, Julius Caesar Scaliger. The

calendar in which 4713 B.C. was speciÖed was the calendar in force in his day,

i.e. the Julian calendar (named after the Roman emperor). In our modern

Gregorian calendar, the same date would be displaced by Öve weeks in time.

Since days are numbered consecutively, if the Julian dates of any two days

are known, a simple subtraction su¢ ces to determine how many days separate

them. Julian dates for a wide range of dates are published, e.g. in the American

22
Ephemeris. Thus, for example, Dec. 21, 1960 is Julian day 2,437,290 while Jan.

1, 1900 is Julian day 2,415,021 and there are 22,269 days between them.

To take advantage of this simplicity in automatic data processing, routines

are required to translate from calendar date to Julian date and back again.

Many such routines have been written and are available from various sources.

The ones we have provided have several features which we feel make them more

than usually useful.

We have extended the concept from days to hours. Thus, a meteorologist

who must think in 3, 6, and 12-hour increments does not have to separately keep

track of days and hours. Furthermore, changing from one time zone to another

can be accomplished by simply adding or subtracting the appropriate number

of hours without concern about possible changes of date. We have deÖned the

Julian Hour to be the number of hours elapsed since 0000 hours, Jan. 1, 4713

B.C. (Julian calendar)2

Another feature is that out-of-range arguments are allowed in the input for

calendar dates, and treated in a logical fashion. This has several beneÖts for

users of the routines, as will be elaborated in the examples below.

The FORTRAN integer function JULHR will calculate the Julian hour for

given information in integer form on year, month, day, and hour. To use, the

programmer should insert the statement

JHR=JULHR(KYR, KMO, KDA, KHR)

in his program. The integer values KYR, KMO, KDA, and KHR which must

be supplied by the programmer represent the year, month, day, and hour, respectively of the calendar date of interest.

Usage notes for JULHR:

1. The arguments are required in order of increasing resolution; i.e. the

largest time unit Örst, the smallest last

2. The year input may be the actual year (e.g. 1980), or it may be a value

from 0 to 99, in which case 1900 will be added before processing (e.g. 78

will be treated as meaning 1978).

3. The range of KMO will normally be from 1 to 12. If a number outside

that range is supplied, a su¢ cient multiple of 12 will be subtracted (or

added) from it and the appropriate number will be added (or subtracted)

to the KYR term. Thus, JULHR(1966,13,5,0) will have the same result

as JULHR(1967,1,5,0).

4. The range for KHR will normally be from 0 to 23, while that for KDA will

normally be 1 through 31. If KHR is out of range, su¢ cient days will be

added to KDA and multiples of 24 subtracted from KHR to Öt. If KDA is

2 By international agreement among astronomers, the Julian day starts at noon, 12 hours

after the corresponding civil date commences. In this way, astronomers do not have to start

a new date in the middle of their working night. Of course, users of these routines are free to

deÖne their starting time as midnight GMT, midnight LST, or whatever is convenient.

23
out of range for any month, the remainder will extend into the following

month. Thus, March 35 will become April 4.

5. Because the values of JULHR in the 20th century are so large (up to

58 106), it requires a computer with a word size of no less than 27 bits

(8 digits) for integer arithmetic. If REAL arithmetic has less than 8 digits

precision (such as single precision on the IBM 360/370) an inadvertent

conversion to REAL will cause serious round-o§ problems. For such a

computer, avoid having JHR appear in an expression with variables of

TYPE REAL or constants in which decimal points or exponential terms

appear.

The FORTRAN callable subroutine DATEX will calculate, from any input

Julian hour, the year, month, day, and hour as integers, as well as the name of

the month as a three character alphameric string (from íJANí to íDECí) and

the name of the day of the week (from íSUNíto íSATí). To use, the programmer

inserts the statement

CALL DATEX(JHR, LYR, LMO, LDA, LHR,

C NAMMON, NAMDAY)

in his program. The integer value JHR (supplied by the programmer) is the

Julian hour of interest; the integer values LYR, LMO, LDA, and LHR (returned by the subroutine) give the year, month, day and hour (in integer form)

corresponding to it. The variables NAMMON and NAMDAY returned by the

subroutine are 3-character left-justiÖed Hollerith strings representing the name

of the month and the name of the weekday, respectively; they may be printed

under an A3 FORMAT speciÖcation.

Usage notes for DATEX:

1. The output arguments are supplied in order of increasing resolution; i.e.

the largest time unit Örst, the smallest last.

2. LYR will be the actual year. Thus, dates in the 20th century will have

LYR in the range 1900 through 1999, not 0 through 99.

3. The range of LMO, LDA, and LHR will be 1-12, 1-31, and 0-23, respectively.

4. NAMMON and NAMDAY may be printed out using an A3 FORMAT.

5. Leap years and leap days are handled according to the rules of the present

day (Gregorian) calendar. Thus, if Feb. 29 is input to JULHR and retrieved through DATEX, the result will be Mar. 1 for the years 1975, 1900

and 1978; Feb. 29 for the years 1976, 1980, and 2000.

Following are some examples of the use of the calendar routines:

24
1. Calculate the Julian day of the year; i.e. the number of days since the

beginning of the year (Jan. 1 = 1 and Dec. 31 = either 365 or 366):

JDAY=(JULHR( KYR, KMO, KDA, KHR) -

C JULHR(KYR, 1, 0, 0)) / 24

(Jan. 0 is actually Dec. 31 of the preceding year).

2. Given the year KYR and the Julian day of the year JDAY, calculate the

calendar date:

JHR0 = JULHR( KYR, 1, 0, 0)

CALL DATEX( JHR0 + 24*JDAY, LYR, LMO, LDA, LHR,

C NAMMON, NAMDAY)

3. Calculate the Julian date of a particular calendar date:

JULDA=JULHR(KYR,KMO,KDA,KHR) / 24

The result should match that in the Ephemeris.

4. As between two given calendar dates, determine which is earlier, and by

how many hours:

KGAP = JULHR(KYR2, KMO2, KDA2, KHR2) -

C JULHR(KYR1, KMO1, KDA1, KHR1)

[IF (KGAP) 10, 20, 30

10 code for case date 2 early

20 code for case dates equal

30 code for case date 2 late

5. Compute mean monthly values on daily data:

DIMENSION VALUE(366)

...

code to evaluate KYR and VALUE

...

KHR0=JULHR(KYR,1,0,0)

DO 50 KMO=1,12

KDA1=(JULHR(KYR,KMO,1,0)-KHR0)/24

KDA2=(JULHR(KYR,KMO+1,0,0)-KHR0)/24

NDAYS=KDA2 - KDA1 + 1

AMEAN=0.

DO 30 KDAY=KDA1,KDA2

AMEAN=AMEAN+VALUE(KDAY)

30 CONTINUE

AMEAN=AMEAN/NDAYS

CALL DATEX(KHR0+24*KDA1, LYR, LMO, LDA, LHR, NAMMON, NAMDAY)

PRINT 40,NAMMON,KYR,AMEAN

40 FORMAT(í MEAN FOR MONTH í,A3,I5,í = í,G10.4)

50 CONTINUE

25
Explanation: KHR0 is the Julian hour beginning the last day of the preceding year (Jan. 0 of this year). KDA1 is the Julian day of the year beginning

the month KMO. KDA2 is the Julian day of the year ending the month KMO

(the zeroth day of month KMO+1). For KMO=12, KDA2 is the last day of

December for the current year, or the zeroth day of the 13th month of this year,

or the zeroth day of the Örst month of the following year. This is a sample of the

programming possibilities inherent in the handling of out of range arguments

for JULHR.

26