A Mathematical Model for the Length of Daylight

[Click here for a sample write up of this project.]

                         ``... it was the season of Light, 
                               it was the season of Darkness ...''
                                             A Tale of Two Cities
                                             C. Dickens

This project will study a mathematical model for the length of daylight as a function of the day of the year and of the latitude. In any place with seasons, the number of daylight hours changes each day, due to the fact that the earth's axis is tilted at an angle of 23 degrees 27'. This tilt causes the north pole, for example, to sometimes face towards the sun (during the Northern Hemisphere's summer, when there are 24 hours of daylight at the north pole) and sometimes away from the sun (during the Northern Hemisphere's winter, when there are 24 hours of night at the north pole).

Like all mathematical models, the one we will talk about involves certain simplifications. It has been derived by a geometrical analysis, in which only the tilt of the earth's axis is taken into account. Other factors affecting the actual number of daylight hours have been ignored. These include the earth's motion around the sun during the course of each day, the change in the earth's orbital velocity as the earth's distance from the sun changes (because the earth's orbit is elliptical), atmospheric refraction, precession of the equinoxes, nutation of the earth's polar axis, the nonzero width of the sun, elevation and other factors which your professors have not been clever enough to think of, but are smart enough to ignore.

Our model is available to you in the function DAYLIGHT, which can be downloaded to your calculator by your TA. [It is documented at the end of this file.] This function gives the number of daylight hours for a location at a given latititude LAT as a function of the time x in days since the beginning of the year. (Before attempting to use DAYLIGHT, you will have to assign a value between 0 and 90 to the variable LAT.)

[1] From personal observation, find the length of the day tomorrow here in Lincoln. (Late risers can consult a newspaper, although they should be aware that newspapers obtain the times of sunrise and sunset from tables prepared by the U.S. Naval Observatory, which in turn are determined from a mathematical model, and may be in error by up to a couple of minutes. In any case, indicate how you obtain your value.) How does this compare to the value given by DAYLIGHT? What is the relative error compared of DAYLIGHT's value (versus the actual value)? (Indicate what data you use as inputs to DAYLIGHT.)

[2] Using your own resources, find the daily rate of change of the length of daylight at this time of year in Lincoln. How did you determine this number? What are its units? [Think carefully about what units you should use in order to make your answer as clear as possible.] What value do you obtain from DAYLIGHT? Explain carefully what calculations you did to determine this value. What is the relative error of DAYLIGHT's value?

[3] Which day or days does DAYLIGHT say are the longest in Lincoln? Which does it say are the shortest? How long does it say they are? Would your answers change if you lived in Miami, FL or Nome, AK?

[4] Make a table of the rate of change of DAYLIGHT as a function of the day. (Be sure to explain what calculations you make.) Also draw a graph of this data against time. On which day(s) is this rate of change the greatest? On which day(s) is it the least? On which day(s) is the rate of change zero? Discuss the connection between these results and those of [3].

[5] During which season(s) of the year can the length of the day in Lincoln be 10 hours? On which day(s) of the year is the day 10 hours long in Lincoln? Suppose the next day there are between 2 and 3 more minutes of daylight; what season must it be? What day must it be? Why?

[6] Can you use a stopwatch to find the latitude? For example, if you measure the number of daylight hours somewhere to be about 16 hours long, where on earth might you be? If the next day is about 5 minutes longer, what day is it and what must your latitude be? Is this answerable? Is there more than one answer? How many are there? Why? How would your being at the equator or above the Arctic Circle affect your ability to answer such questions?

Notes & Advice: There are actually two versions of our model. The function FASTDAY is faster but only works for latitudes below the Arctic circle.

When you present your answers you should write a coherent report that is comprehensible to someone who doesn't have the list of questions in front of them. You should write clearly, in English, using mathematical notation where needed. When you do some calculations make sure you say what you're doing. On the other hand, if you're doing similar calculations repeatedly then you needn't explain the later ones in excessive detail. The sections of the project indicated above needn't be the same as in your presentation. Provided you discuss all the issues raised above, your presentation's organization is up to you.

For Archival Reference Only

The model used in this project regards the number N(t) of hours of daylight as a continuous periodic function with period 365 of a real variable t. It ignores the rotation of the globe; instead it converts t to an orbital position and determines the sunlit fraction f of the circle C of latitude LAT degrees when the Earth is at the specified orbital position. Then N(t)=24*f is the model's output for the inputs t and LAT. (In the project handout above x is used in place of t, but it is convenient in the present discussion to use x for something else.) Comment: Actually, to be more precise, f is the fraction of C in view of the center of the face of the Sun, ignoring atmospheric refraction.

To determine f, consider a (noninertial) coordinate system in which the origin is at the center of the circle C, and the x and y axes lie in the plane P containing C. The z-axis lies on the Earth's polar axis, pointing northward. To fix the x-axis, let l be the line through the center of the Earth, passing through the Northern hemisphere, and parallel to the axis of the Earth's orbit around the sun. Let the x-axis be the line in P through the center of C and through the point p where l and P meet. Choose units such that C has radius 1 and such that the x-coordinate T of p (note that p is on the x-axis) is positive. Then one calculates that T=(tan(23.45*PI/180))*(tan(LAT*PI/180)), where 23.45 is the tilt of the Earth's axis in degrees. The y-axis is now chosen so that the x, y and z-axes comprise a righthanded coordinate system.

Let Q be the plane separating the lit and dark sides of the Earth; thus Q moves throughout the year. Let q be the line where P and Q meet. Then in the x, y, z-coordinate system fixed above, q will be seen to rotate with center p, making one full revolution per year. On the Spring Equinox (i.e., t=82), q coincides with the x-axis; thus the angular measure of the angle (call it a) formed by q and the x-axis is, on any day t, (t-82)*2*PI/365 radians, and we want to determine the fraction f of C which, in the figure below, is above the line q. From the figure and geometry one sees f=(PI+2b)/(2*PI). From the law of sines we have 1*sin(b)=T*sin(a). Thus f=(PI+2*arcsin(T*sin(a)))/(2*PI), with a=(t-82)*2*PI/365 and T as above.

Diagram goes here.

The FASTDAY function referred to above is precisely N(t) as defined above. However, when LAT is above the Arctic Circle, p is outside C and |T*sin(a)| can be bigger than 1, hence outside the domain of definition of arcsin, with the result N(t) does not graph properly on the TI-85 calculator for high latitudes. Thus DAYLIGHT uses conditionals essentially to extend the definition of arcsin so that in addition arcsin(u) is -PI/2 for u < -1 and PI/2 for u> 1. However, these conditionals greatly slow execution on the calculator.

Additional comments:
Note that the angle a used to parametrize the orbit of the Earth around the Sun lies in the plane C. It is more natural but a bit more complicated to use the projection of angle a in the plane of the ecliptic. Call this projected angle A. If we use A to parametrize the Earth in its orbit around the sun we get a new formula: f=(PI+2*arcsin(T*sin(A/(1+cos^2(A)tan^2(23.45*PI/180))^(1/2))))/(2*PI). This is the same as found in the Mathematical Association of America TEAM learning module "A path to applied mathematics: Hours of Daylight", by Prof. James Choike, of Oklahoma State University.

Here are some additional resources related to this topic:
There is the "Almanac for computers 1991 (final edition) ", (in our library its call # is D213.14:991), put out by the USNO Nautical Almanac Office. It gives formulas for hours of daylight accurate for typical latitudes to within a minute or so. For polar latitudes in the periods where the amount of daylight can rapidly change from 24 hours to 0 and vice versa, even this formula will be inaccurate, since it doesn't totally account for the motion of the Earth in its orbit between sunrise and sunset. (However, it would not be hard to fix this by using an iteration based on the basic formula.) In this case the formula for f is: f=[arccos(((sin d) - sin L cos g)/(cos L sin g))]/Pi, where here L is the latitude in radians, d is 50' of arc in radians, and cos g = sin A sin i, where i is 23.45 degrees converted to radians, and A parametrizes Earth's orbit (with A=0 being the vernal equinox, A=Pi/2 being the summer solstice, etc.). To get a good parametrization in terms of time t in years, let s=2tPi/365.25625 - 3.289Pi/180, and substitute A = s + e sin(s) + 1.25e^2 sin(2s) - 77.366Pi/180, where e is the eccentricity of the Earth's orbit (e is about 0.01671). This formula with this parametrization gives pretty good agreement with USNO published values (say within a minute or two at typical latitudes). It still ignores the fact that the Earth moves in its orbit over the course of a day. To get better agreement, you'll want to take this into account, as the formulas in the "Almanac for computers 1991" do (partially) by computing times for sunrise and sunset independently.

For explicit computer code (FORTRAN or BASIC) see USNO Circular no. 171, 2-19-1987, "Computer Programs for Sun and Moon Illuminance with Contingent Tables and Diagrams" (call # D213.4:171).

Other resources are also available at the USNO website. There is available (for sale) MICA (Multiyear Interactive Computer Almanac) at http://riemann.usno.navy.mil/AA/ Also visit http://www.whnt19.com/weather/sun_set.htm for on-line calculations.