Note: A bug in the conversion program has subsituted the greek letter "rho" for vector arrow signs

Using a Digital Terrain Model to Calculate Visual Sunrise and Sunset Times

Chaim Keller

Midrash Bikurei Yosef

P.O. Box 35078, Jerusalem 91350, Israel

FAX: 972-2- 5713765

John K. Hall

Marine Geology Division, Geological Survey of Israel

30 Malkhei Israel Street, Jerusalem 95501, Israel

e-mail: johnhall@mail.gsi.gov.il

Abstract--Since its advent the Digital Terrain Model (DTM) has been employed widely in the sciences for the solution of problems requiring a digital model of landforms. In this paper we describe a new use of the DTM in the calculation of highly accurate visual sunrise and sunset times that are required by the observant population in Israel. We have employed ray tracing to determine the effect of atmospheric refraction through a simplified layered atmosphere. A general analytic expression for the atmospheric refraction was determined from these calculations as a function of the observer’s height for two model atmospheres known as the subtropical summer and winter atmospheres . These expressions determine the general magnitude of the refraction as a function of apparent view angle of the observer. We also determined a simplified analytic expression for the effect of atmospheric refraction on the vertical angular profile of the mountainous horizon as calculated from the DTM (atmospheric refraction magnifies the mountainous features of the horizon). These expression are then used in calculating the apparent vertical angular position of the sun as a function of time. The time when the upper limb of the sun first (last) appears to rise (set) over the horizon adjusted for the effect of refraction determines the time of the visible sunrise (sunset). Comparison with observations have shown that the visible sunrise and sunset times can be typically determined to better than 15 seconds using the 25 meter DTM of Israel and eastern Jordan.

Key Words: DTM. DEM. Refraction. Sunrise. Calendar.

The visible sunrise and sunset times are defined as those times at which the sun rises and sets, respectively, over the mountainous horizon. This changes daily as the sun rises and sets over different terrain as it progresses in its apparent daily motion. The calculation of the visible sunrise and sunset times is therefore impossible without an accurate topographical model of the eastern and western horizons, respectively. In the past this difficulty could be overcome only through daily observations or by optical measurements with a theodolite. Since both of these methods are laborious they can only be applied in special situations. On the other hand, the digital terrain model (DTM) provides a flexible and rapid means for obtaining an analytical description of the visible horizon. This in turn can be used for computing the sun’s visible sunrise and sunset times.

Neglecting atmospheric refraction, if h₁ and h₂ are the heights of the observer and of the obstructions at a distance, , from the observer (see Figure 2), then the obstruction will appear to make an angle, , with the horizontal at an azimuth angle b, where

, (1)

(2)

and,

, h_i, and (3)

and,

, , (4)

; , (5)

where , are the latitude and longitude respectively (i = 1 for the observer, and i = 2 for the obstruction in the distance). The physical shape of the eastern and western horizons are determined by finding the maximum viewing angle, , for any particular azimuth, b, among all the relevant elevation points, (, h_i ), contained in the DTM of Israel and eastern Jordan (Hall, 1996; 1997b). This DTM, depicted in Figure 1, consists of elevations on a 25 meter grid, derived from the 1:50,000 scale mapping using methods described in Hall et. al. (1990).

The visible horizon is not determined solely by physical landforms. Rather, atmospheric refraction is an important contributor to the shape of the observable horizon if the distances to the landforms that determine the horizon are significantly greater than 15 km (as shown by our calculations). If is the viewing angle of the mountainous horizon after atmospheric refraction then

, (6)

where is given by Eq. 1, and is a positive correction factor. depicts the magnifying effect of atmospheric refraction on physical landforms on the horizon. More generally, the apparent viewing angle of the sun’s center, , is related to the total atmospheric refraction, , as follows,

, (7)

where is the angle elevation of the sun’s center with respect to the horizon without atmospheric refraction. can be calculated to a sufficient approximation by using highly idealized layered model atmospheres such as employed by Menat (1980). (Unless mentioned otherwise, the optical-atmospheric physics in this work is based on the work of Menat (1980) and the literature cited therein). The temperature and pressure dependencies of the layered atmospheres are given by

=+ h; (8)

= (9)

The values of , , , and are given in Table I for the two model atmospheres used in this work. They are known as the subtropical summer and winter atmospheres (Selby, 1976, Valley, 1965) in accordance with their characteristic temperature and pressure dependencies with altitude. The atmospheric refraction, , can be calculated for any desired set of light trajectories corresponding to observation angles, , heights h₁, h₂, and distances D = (see Figure 2). This is done by solving Snell’s law for arbitrarily small steps along the path of the light trajectory that corresponds to the above boundary conditions, i.e.,

(10)

where is the angle the light trajectory makes to the normal of the atmospheric layer i. Neglecting the effect of humidity, the index of refraction, , at any altitude, , is given by

/ l² (11)

; ; l = 0.6 mm

where the height dependencies of the pressure, , and temperature, , are given in equations 8, 9 and in Table I. In practice we calculated the true viewing angle at the observer, , and the total atmospheric refraction, , separately. We performed these apparently redundant calculations in order to determine an analytic expression for the correction factor (Eq. 6) to the view angle, (Eq. 1), as describebelow.

The total atmospheric refraction, , was calculated by ray tracing for a range of observer heights, h₁, and view angles, . In calculating , we traced the path of the light in the model atmospheres by solving equations 10-11 in 2 meter steps from the observer to the outer limit of the model atmospheres. [For negative view angles, , the path of the light ray was calculated in steps of minus 2 meters until the light ray leveled off locally (a = 90 degrees) due to the curvature of the earth. After this point the path of the light was calculated in +2 meters steps until the light reached the outer limit of the model atmosphere. is the sum total of the refraction from these two regions.] The results for could be adequately fitted to a second order polynomial in (degrees) that is also a function of the observer height, h_1, viz.,

(degrees). (12)

The values of the polynomial coefficients for the observer heights, h₁ (with respect to sea level), relevant to this work are given in Tables II-III. Figure 5 plots the values of for 750 h₁ 850 (meters) for the winter atmosphere.

In order to calculate (see Eq.. (6) ) the path of the light rays was traced backward from the observer to the obstruction. A search routine was written to determine at which angle the ray just glanced off the obstruction at distance D at height h₂. The search began at angles , and stepped in angle until the ray just glanced the obstruction at distance D at height h₂. The angle, , which meets this condition is the solution to Eq. 6 for the chosen value of h₁, h₂, and D. The results of these calculations were fitted to a simplified analytic expression given in Eq. 13. Due to the small difference in between the summer and winter atmospheres, it was sufficient to use the average , , for either atmosphere.

(13)

where D is in km, and h₁, h₂ are in meters, is in degrees, and

.

The small discontinuity in near in Eq. 13 and in Figure 3 is caused by using a simplified analytic expression for the refraction that models the behavior of as two intersecting straight lines. The simplified expression gives values of that differ no more than by 0.02^o from the actually calculated values. This produces a negligible error of 6 seconds in sunset (sunrise) times for our latitude. The results are plotted in Figures 3-4 for some typical values of D, h₁, and h_2.. Note that in Figure 3 the refraction correction, , is asymmetric after the crossing point at . This is probably due to light bending more in order to glance off the obstruction as h₁ becomes larger than h₂ (see Figure 2). As seen in Figure 4, the further the obstruction, the larger the refraction correction, , becomes. This is due to the light trajectory bending more for larger distances, D, while for very short distances ( D 15 km ) the trajectory is constrained to be a straight line and = 0, as can be seen in Figure 2.

The visible sunrise and sunset times were calculated by determining at which time the upper limb of the sun passed over the mountainous horizon. A search program was written for this purpose. This program begins looking for the sunrise or sunset time over the mountains at the time of the sunrise or sunset for an observer at height h₁ without mountains in the distance (known as the astronomical sunrise or sunset times). The program then advances, or goes backward, in 2.5 second steps in time until the visible sunrise or sunset times are reached, respectively. If _sun, is the apparent angular altitude of the upper limb of the sun and a_alt is the angular altitude of the center of the sun without refraction (whose value is determined by the equation of time, see, e.g., Astronomical Almanac, USGPO, 1996), then

_sun = a_alt + e_sun + (14)

where e_sun is the angular half width of the sun (0.2666 degrees) and is the total atmospheric refraction as given in Eq. 12. Since itself depends on _sun, Eq. 14 must be solved by an iterative process. For an initial guess for , we use the total astronomical refraction at the astronomical sunrise or sunset as calculated from ray tracing. We then introduce the result of Eq. 14 into Eq. 12 and reiterate Eq. 14 with the new value of until we obtain a constant value of _sun to the accuracy of the calculations. We have found that Eq. 14 converges rapidly after only 1 to 4 iterations. This process is repeated for each step in time until the apparent altitude of the upper limb of the sun, _sun, at the current solar azimuth is equal to (Eq. 6) at that same azimuth. At this point the visible sunrise or sunset time for that day has been reached.

Figure 6 shows the result of our calculation (black squares) of the visible sunrise times for the Armon Hanatziv neighborhood of Jerusalem (35^o 14’ 18.5” longitude east, 31^o 44’ 59.8” latitude north), h₁ = 754.9 m. The sunrise horizon is defined by the Moav mountains of Jordan 40-70 km to the east whose peak altitudes, h₂, range from 700 to 1200 meters. Ten years of near-daily observations by Rabbi Avraham Druk are superimposed over our calculated times in Figure 5. We have used the summer atmosphere for days 85 to 290 and the winter atmosphere for the rest of the 365 days of the year. The 10 years of observations have been suitably shifted in order to neutralize the effect of time shifts between different years (in most part caused by the 4 year cycle). The results show that the sunrise time can be calculated to an accuracy of about 15 seconds for more than 90 percent of sunrise times observed. Comparison of observations to our calculations for other places in Israel have been consistent to the same accuracy. It is our belief that adding a term for humidity in the index of refraction (Eq. 11) can improve the accuracy of our calculations for the winter months. Figure 7 gives an example of the sunrise calendars that we provide through a country-wide publication (The Bikurei Yosef Tables of the Visible Sunrise and Sunset Times).

For a distant horizon, as that of Figure 6, elevation inaccuracies in the DTM ( meters maximum, but typically no worse than meters) produce a quite negligible effect on the calculated sunrise times. However, these inaccuracies will make the major contribution to the uncertainty in the sunrise and sunrise times for distances less than 5 km. We therefore cannot provide tables for those cities that have mountainous obstructions within 5 km.

The observations shown in Figure 6 are in contradiction to the cof Schaefer and Liller (1990) who cit is impossible to the sunrise times to better than 4 min. They obtained this conclusion by averaging observations from varying heights. However this neglects the dependence of the rms deviation of the sunrise times on the length of the light’s air path. For example, the length of the light’s air path to the horizon, L, for an observer of height, h₁, is (Menat, 1980, his Eq. 41),

L . (15)

It can be shown that a proper analysis of their data, taking account of Eq. 15, is consistent with the observations of sunrise times for Israel.

Acknowledgments--The authors would like to acknowledge their deep gratitude to Dr. Dror Bendavid of Midrash Bikurei Yosef for initiating and supporting this project. We would also like to thank Dr. Mordechai Menat for providing us with his computer codes and for his useful suggestions and comments. Finally we would like to thank Rabbi Avraham Druk for kindly providing us with his sunrise observations. The 25 meter DTM of Israel is the copyrighted property of the Survey of Israel.

Astronomical Almanac, 1996. U.S. Government Printing Office, p. C24.

Hall, J. K., 1996. Topography and Bathymetry of the Dead Sea Depression, Tectonophysics, 266 (1-4), pp. 177-185. This is one of thirteen papers in the Editor’s Personal Selection of the top papers published in Tectonophysics over the past two years. It appears in full on their website at http://www.elsevier.com/locate/tecto.

Hall, J. K., 1997a. Three-dimensional map of Israel and vicinity. Grayscale sheet of land topography at 1:500,000 scale with physiographic provinces by E. Zilberman and R. Bogoch (English Edition). Geological Survey of Israel, Jerusalem, April 1997.

Hall, J. K., 1997b. Topography and bathymetry of the Dead Sea Depression. In: The Dead Sea - The Lake and its Setting, eds. Tina M. Niemi, Zvi Ben-Avraham, and Joel R. Gat, Oxford Monograph on Geology and Geophysics No. 36, Oxford University Press, New York, pp 11-21.

Hall, J. K., Schwartz, E., and Cleave, R. L. W., 1990. The Israeli DTM (Digital Terrain Map) Project. In: Microcomputer Applications in Geology, Volume II, eds. J. Thomas Hanley and Daniel F. Merriam, Pergamon Press, pp. 11-118. Includes FORTRAN-77 DTM program on diskette as part of the COGS (Computer Oriented Geological Society) public-domain library.

Menat, M., 1980, Atmospheric Phenomena before and during Sunset, Applied Optics, 19, 20, pp. 3458-3468.

Selby, J. E. A., et. al., Atmospheric Transmission from 0.25 to 28.5m. Supplement LOWTRAN 3B, In: AFGL Report TR-76-028, Air Force Systems Command, USAF, 1976.

Valley, S. L., Ed., Handbook of Geophysics and Space Environments, McGraw-Hill, New York, 1965, Chaps. 7 and 9.

Figure 1 Shaded relief image of the study area based upon the 25 meter DTM of Israel and eastern Jordan (Hall 1977). Jerusalem is located in the relatively flat area in the mountains north west of the Dead Sea. For the high places of Jerusalem, sunrise appears over the Moav Mountains in eastern Jordan, east of the Dead Sea.

Figure 2 Depicts the path of light, in the model atmospheres, glancing off the obstruction at height, h₂, and finally reaching the observer at height h₁. In the figure, the obstruction appears to make a negative angle, , with the observer’s local horizontal. Without refraction, the obstruction would make an angle, <. Other variables are defined in the text. The dimensions are not to scale.

Figure 3 Plot of the refraction correction, , to the view angle, , as a function of the height of the observer, h₁ for an obstruction height h₂ = 500 m, and a distance to the obstruction, D = 50 km. (See text for details.)

Figure 4 Plot of the refraction correction, , to the view angle, , as a function of the distance, D, to the obstruction, for an observer height h₁=500 m, and an obstruction height, h₂ = 500 m. (See text for details.)

Figure 5 Plot of the total atmospheric refraction, , as a function of view angle, , for the winter atmosphere for 750 m h₁ 850 m.

Figure 6 Calculated sunrise times (black line) superimposed over 10 years of near daily observations of the sunrise by Rabbi Avraham Druk from the Jerusalem neighborhood, Armon Hanatziv. (See text for details.)

Figure 7 Example of the Bikurei Yosef Tables of the Visible Sunrise and Sunset for the city of Jerusalem. The tables are calculated for the city of Jerusalem according to the Hebrew Calendar of the year 5760 which spans the calendar years 1999-2000. The year 5760 is a leap year having 13 months of length 29 or 30 days. The tables are read from top to bottom, right to left. The first time entry (the Jewish Near Year) corresponds to September 11, 1999. As explained in the text, we are over 90% confident that the calculated times are accurate to better than 15 seconds. Jewish religious considerations therefore require us to add 15 seconds to each sunrise time, and subtract 15 seconds from each sunset time. The sunrise and sunset times are then rounded to the nearest upper or lower 5 seconds, respectively.

Figure 2

Figure 3

Figure 4

Figure 5

Table 1 Value of coefficients in Eq. 12

Table 2 Values of coefficients in Eq. 12 for the winter atmosphere

Table 3 Values of coefficients in Eq. 12 for the summer atmosphere

h1	a(h1)	b(h1)	c(h1)	d(h1)	e(h1)	f(h1)
h1 1150 meters	0.16878		0.0148	0.0001	0.47806	0.07447
1050 h1 < 1150	0.16901		0.01575	0.00007	0.4834	0.07641
950 h1 < 1050	0.16896		0.02169	-0.00014	0.52128	0.09025
850 h1 < 950	0.16838		0.02111	-0.00012	0.5161	0.08825
750 h1 < 850	0.16876		0.0232	-0.00019	0.52786	0.09266
650 h1 < 750	0.16907		0.02383	-0.00021	0.5303	0.09358
550 h1 < 650	0.1691		0.02976	-0.00042	0.56622	0.10671
450 h1 < 550	0.1696		0.02905	-0.00039	0.56005	0.10422
350 h1 < 450	0.1691		0.01416	0.00011	0.46823	0.06829
250 h1 < 350	0.169		0.01168	0.0002	0.45241	0.06206
150 h1 < 250	0.16941		-0.01107	0.0083	0.32072	0.0059
50 h1 < 150	0.16924		-0.00922	0.0008	0.33125	0.1087
-50 h1 < 50	0.16872		0.0305	-0.0021	0.55065	0.10908
-150 h1 < -50	0.16872		0.03083	-0.00022	0.55083	0.10908
h1 < -150	0.16870		0.03119	-0.00022	0.55118	0.10915

Table 2