During the course of one year, I photographed 95 sunrises along the east bank of the Columbia River from my home on the west bank in Richland, WA (46.3° N latitude). I then calculated the seasonal phenomena listed in the title with the intention of explaining the variation of the photographed sunrises. The calculations use a simplified model of the Sun-Earth system and employ rotation matrices to predict the path of the Sun, as observed at any location in the northern hemisphere, throughout the year. These predictions are in good agreement with those listed by NOAA and also with the photographic data. The analysis presented here provides a novel way to calculate and understand the seasonal variations of visible sunlight.

## I. INTRODUCTION

Several years ago, I began a retirement project by taking daily photographs of sunrises visible from my home on the Columbia River in Richland, WA, during an entire year. Figure 1 shows a sample of those photographs taken in September and October 2018. The photos and their analysis were published in Ref. 1. More recently,^{2} the changing number of daylight hours and the elevation of the Sun at noon were studied using geometric diagrams.

Other authors treat the same phenomena. The text by Meeus^{3} discusses corrections to the basic theory due to small effects, such as the refraction of sunlight in the Earth's atmosphere, that affect the number of daylight hours. References 4 and 5 are popular texts for amateur astronomers. Levine^{6} used vector analysis to understand the duration of daylight and the position of the rising and setting Sun.

In contrast to the above treatments, this paper studies the apparent motion of the Sun by employing *rotation matrices*. Using a simplified model of the solar system, discussed in Sec. II, this approach sacrifices some accuracy to provide physics students and instructors with an easier analysis and deeper understanding of the phenomena. Reference 2 provides a good background to this approach.

## II. SIMPLIFIED MODEL AND THE OP-COORDINATE SYSTEM

Our simplified model starts with a heliocentric coordinate system: The Sun is located at the origin of the orbital plane (OP) coordinate system, in which a spherical Earth circles the Sun in the *x*_{OP}-*y*_{OP} plane. One hemisphere of the Earth is illuminated by sunlight, represented by a unit vector **S** that is parallel to the line joining the centers of the Sun and Earth; the other hemisphere is in darkness. Most of the calculations in this paper will use the latitude of the author's home, which is 46.3° N.

Rather than using this heliocentric system, it is easier to proceed using the geocentric coordinate system shown in Fig. 2 with the Earth at its origin and the Sun revolving around it in a circular path in the OP plane. The vector **S**, when directed toward the Earth's center, makes an angle $\kappa $ with the negative *x*_{OP} axis, where $\kappa $ ranges from 0° to 360° during one year. The asterisks in Fig. 2 indicate the circular boundary between the light and dark hemispheres, called the light-dark circle (LDC). As **S** changes the direction each day, the plane of the LDC changes its orientation, causing the number of daylight hours at a particular location on Earth to change as well.

*κ*= 0° and DN = 1. While the angle

*κ*continually increases, the simplified model assumes that it remains constant during each day and increases abruptly by $ 360 \xb0 / 365 = 0.986 \xb0$ at the beginning of each new day

**S**lies along the

*x*

_{OP}axis, and the LDC lies in the

*y*

_{OP}–

*z*

_{OP}plane. Points on the LDC (asterisks in Fig. 2) are given in terms of OP coordinates by

*η*ranges from 0° to 360° relative to the +

*y*

_{OP}axis.

When $\kappa \u22600$, it is helpful to construct the vector **F** that points from the Earth's center to a specific location on the LDC. It has a horizontal component $ F h= R E \u2009cos\u2009 \eta $ and a vertical component $ F v= R E \u2009sin\u2009 \eta $ that lie in the plane of the LDC. Let us imagine that the LDC is rotated from its initial orientation $\kappa =0$ to the orientation shown in Fig. 2. In that case, the vertical component $ F v$ remains unchanged but the horizontal component $ F h$ now has components along the –*x*_{OP} axis and the *y*_{OP} axis. The equations for the LDC in the OP-coordinate system are now given by

## III. THE E-COORDINATE SYSTEM

### A. Equation of the LDC in E-coordinates

We now introduce the E-coordinate system shown in Fig. 3, where the *z*_{E} and *x*_{E} axes are rotated by 23.5° relative to the OP-coordinate axes, and *z*_{E} is aligned with the Earth's N-S axis. The *x*_{E} and *y*_{E} axes lie in the planet's equatorial plane. Note that the Earth rotates with respect to this coordinate system; we will later introduce a coordinate system (called Q) that is fixed with respect to points on the Earth's surface.

The locations along the LDC that are visible in this figure are all experiencing sunset at the moment shown, since an observer rotating counter-clockwise with the Earth is travelling from sunlight into darkness. On the opposite side of the planet, hidden from view in Fig. 3, locations along the LDC are experiencing sunrise since they are traveling from darkness into sunlight. To find the number of hours of daylight and to understand the apparent shift in location and direction of sunrise shown in Fig. 1, we need to determine the *x*_{E} and *y*_{E} coordinates of the sunrise and sunset points for any day of the year and for any latitude *L*.

*θ*= 23.5° about the $ y OP$ axis, Eq. (A4) implies

*x*(Eq. (A6)) and

*z*axes (Eq. (A8)).

### B. Transformation of S from the OP- to the E-coordinate system

**S**$$ always lies in the

*x*

_{OP}–

*y*

_{OP}plane and has OP coordinates given by $ cos \kappa , sin \kappa , 0$. Using Eq. (6) with $\theta =23.5\xb0$ (the Earth's tilt relative to its orbital plane), the components of

**S**can be transformed from the OP- to the E-coordinate system, with the following result:

*S*x

_{E}and

*S*y

_{E},

**S**xy is the component of **S** in the *x*_{E}–*y*_{E} plane. In a three-dimensional diagram, vector **S** is parallel to the line joining the center of the Sun to the center of the planet.

## IV. SUNRISE AND SUNSET IN E-COORDINATES

In Fig. 4, the 46.3° N latitude circle and the LDC are plotted in E-coordinates for February 20 (DN = 62, *κ* = 60°). The intersections between the LDC and the latitude circle correspond to sunrise and sunset locations. For *κ* = 60°, the OP coordinates of the LDC are found from Eqs. (5): $ x O P , \u2009 y OP , \u2009 z OP \u2009= R E( \u2212 3 / 2 \u2009 cos \u2009 \eta , \u2009 1 2 cos \u2009 \eta , sin \u2009 \eta )$. These are then substituted into Eq. (6) to find their E components

Intersections with the 46.3° N latitude circle are easily found, because for that circle *z _{E}*/

*R*= sin(46.3°) = 0.7230. A simple way to proceed is to search for solutions to Eqs. (9) using a spreadsheet. For $\kappa =60\xb0$, the results are shown in Table I. Clearly, $ z E / R E\u22450.7230$ when

_{E}*η*≅ 68° and 153°. If more accuracy is desired, these two angles can be used as initial estimates to locate the roots of Eq. (9c) numerically, with the results

*η*= 68.2° and 153.1°. Evaluating Eqs. (9a) and (9b), the E-components are found to be (–0.6652, 0.1857, 0.7230)

*R*and (0.5278,–0.4463, 0.7230)

_{E}*R*for sunrise and sunset, respectively. The sunrise and sunset locations for any other latitude

_{E}*L*and angle

*κ*can be found in the same way. A second method using algebra to find the E coordinates of sunrise and sunset is described in the supplementary material.

^{9}

LDC . | Point on LDC in OP coordinates . | Point on LDC in E coordinates . | ||||
---|---|---|---|---|---|---|

η (°)
. | x_{OP}
. | y_{OP}
. | z_{OP}
. | x_{E}
. | y_{E}
. | z_{E}
. |

67 | –0.3384 | 0.1954 | 0.9205 | –0.6774 | 0.1954 | 0.7092 |

68 | –0.3244 | 0.1873 | 0.9272 | –0.6672 | 0.1873 | 0.7209 |

69 | –0.3104 | 0.1792 | 0.9336 | –0.6569 | 0.1792 | 0.7324 |

⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ |

152 | 0.7647 | –0.4415 | 0.4695 | 0.5140 | –0.4415 | 0.7354 |

153 | 0.7716 | –0.4455 | 0.4540 | 0.5266 | –0.4455 | 0.7240 |

154 | 0.7784 | –0.4494 | 0.4384 | 0.5390 | –0.4494 | 0.7124 |

LDC . | Point on LDC in OP coordinates . | Point on LDC in E coordinates . | ||||
---|---|---|---|---|---|---|

η (°)
. | x_{OP}
. | y_{OP}
. | z_{OP}
. | x_{E}
. | y_{E}
. | z_{E}
. |

67 | –0.3384 | 0.1954 | 0.9205 | –0.6774 | 0.1954 | 0.7092 |

68 | –0.3244 | 0.1873 | 0.9272 | –0.6672 | 0.1873 | 0.7209 |

69 | –0.3104 | 0.1792 | 0.9336 | –0.6569 | 0.1792 | 0.7324 |

⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ |

152 | 0.7647 | –0.4415 | 0.4695 | 0.5140 | –0.4415 | 0.7354 |

153 | 0.7716 | –0.4455 | 0.4540 | 0.5266 | –0.4455 | 0.7240 |

154 | 0.7784 | –0.4494 | 0.4384 | 0.5390 | –0.4494 | 0.7124 |

*R*. Note that $ S S \u2192= S R \u2192=0.6909 R E= R E cos \u2009 46.3 \xb0 .$ The angle

_{E}*β*between these two vectors is found using the dot product

Figures 5(a)–5(h) show the distribution of sunlight and darkness in the northern hemisphere for eight dates spaced throughout the year. In each figure, the $ x E , \u2009 y E$ coordinates of the LDC (asterisks) are plotted for $0\u2264\eta \u2264180\xb0$, i.e., points in the northern hemisphere, using Eqs. (6). In the upper left corner for each figure, the angle *κ* is given, followed by the date and the day number (DN). In the upper right corner, the angle *φ* (defined above) and the number of daylight hours are shown. The vector **S**_{xy}, defined in Eq. (8), is drawn for each figure as the perpendicular bisector of the line connecting the sunrise and sunset points.

Many seasonal features can be understood using Fig. 5, such as: (1) the equator always has 12 h of daylight and 12 h of darkness; (2) for the Arctic Circle ( $L=90\xb0\u201323.5\xb0=66.5\xb0$), there is total darkness on December 21 and total sunlight on June 21; (3) there are 12 h of daylight everywhere on Earth on the equinoxes (March 21 and September 21); (4) in the northern hemisphere, the hours of daylight increase from December 21 (winter solstice) to June 21 (summer solstice), and then decrease until December 21; (5) for any date, as you travel to lower latitudes, the hours of daylight increase. Notice (6) also that from December 21 to January 19, the hours of daylight increase only 0.6 h, while from February 20 to March 21 they increase by 1.6 h. Thus, to an observer on Earth, the Sun appears to travel slowly at the solstices and faster at the equinoxes.

Figure 5 of Ref. 2 displays the number of hours of daylight, calculated as above, vs. the date (or day number), and compares them to information obtained from the National Oceanic and Atmospheric Administration (NOAA) website.^{10} The precise NOAA values are about 0.2 h longer than those calculated here, a difference of 2.4% on December 21 and 1.3% on June 21. The discrepancy is due to small effects that are not considered here, such as the refraction of sunlight in the Earth's atmosphere.

## V. THE OBSERVER AND THE Q-COORDINATE SYSTEM

In what direction does an observer on Earth look to view the Sun? To answer this question, we introduce a new coordinate system Q that rotates with the observer about the planet's N–S axis. As shown in Fig. 6, the + $ y Q$ axis points north, the $ + x Q$ axis points east, while the $ + z Q$ axis is in the observer's vertical direction. At noon on December 21, the sunlight vector **S** lies in the plane of the paper, and a straightforward analysis using plane geometry is possible. For other locations of the observer or other times and days of the year, rotation matrices show what an observer sees while watching the path of the Sun. However, keep Fig. 6 in mind, as it illustrates the essential mathematics for determining the components of **S** in the observer's reference frame.

*x*

_{E}axis, so

*κ*= 0. The observer views the Sun in the direction of

**SO**

_{Q}, which is antiparallel to

**S**, i.e.,

**SO**

_{Q}= –

**S**. The dashed line in Fig. 6 is the (local) horizontal axis –

*y*

_{Q}, which points due south. In the figure, we see that the sum 46.3° + 23.5° + E = 90°, so the elevation angle E of the Sun is equal to 20.2° and

**SO**

_{Q}is given by

In Fig. 6, the *observer's* Q-coordinate system (called the “upper” Q system) has its origin at the location of the observer. However, when using a rotation matrix to transform from one coordinate system to another, the origins of the two systems must coincide. Therefore, we define another Q-coordinate system (called the “lower” Q system) whose axes are parallel to the upper Q system but whose origin is located at the Earth's center. The components of **S** are *identical* in the upper and lower Q systems. Similarly, the *x*_{Q} and *y*_{Q} values of the LDC are identical in the two systems, but the *z*_{Q} values differ by the radius of the Earth $ R E$. While we visualize events in the upper Q-coordinate system, we employ the rotation matrices using the lower Q system.

## VI. USING THE Q-COORDINATE SYSTEM TO DESCRIBE THE SUN'S PATH

### A. Vectors S and SO in the Q-system

Two rotations are required to transform from the E-system to the Q-system. Recall that all coordinate systems have their origin at the center of the Earth. The first is shown in Fig. 7(a) for December 21: The –*y*_{E} axis is rotated by angle *θ* about the *z*_{E} axis so that it passes directly *BELOW* the location of the observer. This rotation yields an intermediate coordinate system, which we will call the M-coordinate system and label the rotation matrix **R**_{ME}. To obtain the value of *θ*, we use Fig. 5(a), noting that angle $\beta =125.9\xb0$ and also that the sunlight is symmetrical about the *x*_{E} axis. The positions of sunrise and sunset are then $180\xb0\u2009\u2213 \beta / 2=(117.05\xb0,\u2009242.95\xb0)$. To rotate the −*y*_{E} axis to the sunrise point, $\theta =117.05\u2212270\xb0=\u2212152.95\xb0,$ where the negative angle denotes a CW rotation. The rotation matrix **R**_{ME} is obtained by substituting *θ* = −152.95° in Eq. (A8) of the Appendix.

*z*

_{M}axis, the rotated axis

*y*

_{M}, and the sunrise point. After rotating by $+43.7\xb0$ about the

*x*

_{M}axis, the

*z*

_{M}axis passes through the sunrise point. The rotation matrix

**R**

_{QM}is obtained by substituting +43.7° for angle

*θ*in Eq. (A6) of the Appendix. This yields the Q-coordinate system shown in Fig. 7(b), where the +

*x*

_{Q}(and +

*x*

_{M}) axes are perpendicular to the plane of the paper, pointing outward. The “upper” Q-system can then be constructed at the location of the observer, with the desired arrangement of axes shown in Fig. 6. For a different latitude

*L*, the rotation angle used in

**R**

_{QM}is 90° –

*L*. Of course, the locations of sunrise and sunset would need to be determined to find angle

*θ*. The transformation from E- to Q-coordinates may, thus, be written in matrix notation as $ S Q= R QM\u2009 R ME\u2009 S E$, and the generalized equation for the direction to the Sun, according to the observer, is

*L*is the observer's latitude,

*θ*is the angular rotation of the –

*y*

_{E}axis about the

*z*

_{E}axis, 90 –

*L*is the rotation about the

*x*

_{M}axis, and

*κ*depends on the day of the year.

### B. Sunrise and Noon on December 21

We can now use Eq. (14) to find **SO**_{Q} when the observer views the Sun at sunrise on December 21 (*κ* = 0). Evaluating Eq. (14) for *L =* 46.3° and *θ* = –152.95° yields $ SO Q= 0.8167 , \u2212 0.5770 , \u2009 0.$ The positive sign of SO_{Qx} and the negative sign of SO_{Qy} show that the observer is facing southeast (components to the south and east) to view the Sun. Specifically, the angle *ϕ* between the +*x*_{Q} axis and the vector **SO**_{Q} is given by $ tan \u2212 1 \u2212 0.5770 / 0.8167 = \u2212 35.2 \xb0$. As expected, the *z*_{Q} component is zero at sunrise.

At noon on December 21, the observer is midway between sunrise and sunset, or on the –*x*_{E} axis, as shown in Fig. 7(a). Therefore, the –*y*_{E} axis must be rotated CW by $\theta =\u221290\xb0$ about the *z*_{E} axis. The values for *L*, *κ*, and **S** remain the same as those for sunrise. Evaluating Eq. (14) yields $ SO Q= 0 , \u2212 0.9385 , \u2009 0.3453,$ indicating that the vector points due south. The elevation angle E is the angle between the Sun and the horizon and is given by $E= tan \u2212 1 0.3453 / 0.9385=20.2\xb0.$ This is the value that was found in Fig. 6 using geometry.

### C. Path of the Sun on December 21, March 21, and June 21

#### 1. Path of the Sun on December 21

The apparent path of the Sun on December 21 is obtained by evaluating Eq. (14) throughout the day, as the observer moves from sunrise to sunset. Since $\theta =\u2212152.95\xb0$ at sunrise, and $\theta =\u221227.05\xb0$ at sunset (Fig. 7(a)), the value of *θ* at equally spaced intervals during the day is easily found. A three-dimensional plot of (SO_{Qx}, SO_{Qy}, SO_{Qz}) for 17 times (shown as + signs) between sunrise and sunset is displayed in Fig. 8(a).

#### 2. Path of the Sun on March 21 and June 21

The calculations of **SO**_{Q} for March 21 and June 21 (Figs. 8(b) and 8(c)) are similar to those for December 21, but since the sunrise and sunset locations are different for each date, the value of *θ* in Eq. (14) also changes. In Figs. 8(b) and 8(c), the path of the Sun is shown for these dates. Table II summarizes the calculations for the vector **SO**_{Q} on the three dates shown in Fig. 8. Elevations obtained from the NOAA website^{10} for the three dates are as follows: 20.31°, 44.17°, and 67.15°. The values obtained here are, thus, in excellent agreement with those computed by NOAA.

Date and time . | SO_{Q}
. | Angle ϕ or E
. |
---|---|---|

Dec 21 sunrise | (0.8167, −0.5770, 0) | ϕ = −35.2°, south-east |

Dec 21 noon | (0, −0.9385, 0.3453) | E = 20.2°, due south |

March 21 sunrise | (1.0000, 0, 0) | ϕ = 0, due east |

March 21 noon | (0.−0.7230, 0.6909) | E = 43.7°, due south |

June 21 sunrise | (0.8167, 0.5770, 0) | ϕ = 35.2°, north-east |

June 21 noon | (0, −0.3875, 0.9210) | E = 67.2°, due south |

Date and time . | SO_{Q}
. | Angle ϕ or E
. |
---|---|---|

Dec 21 sunrise | (0.8167, −0.5770, 0) | ϕ = −35.2°, south-east |

Dec 21 noon | (0, −0.9385, 0.3453) | E = 20.2°, due south |

March 21 sunrise | (1.0000, 0, 0) | ϕ = 0, due east |

March 21 noon | (0.−0.7230, 0.6909) | E = 43.7°, due south |

June 21 sunrise | (0.8167, 0.5770, 0) | ϕ = 35.2°, north-east |

June 21 noon | (0, −0.3875, 0.9210) | E = 67.2°, due south |

## VII. COMPARISON BETWEEN CALCULATIONS AND PHOTOGRAPHIC DATA

In this section, we shall compare sunrise photographs taken by the author, several of which are shown in Fig. 1, to the theory developed above. A simplified bird's eye view of the area near the author's home is shown in Fig. 9. Distances along the horizon were obtained by comparing features on Google Maps to the same features visible in the photographs, as discussed in Ref. 1. Figure 9 depicts seven sunrises chosen about one month apart. The asterisks show where the Sun appears at dawn along the bluff of the river. The dotted lines are drawn using the calculated values of *ϕ* for these dates. They converge at point C, which is ≈1320 m west of line AB and is taken to be the origin of the Q-coordinate system. The distance of each sunrise from that on the equinox (March 22) at point D is given by $ D photo=DC\u2009 tan \u2009 \varphi $.

To explore this further, we analyzed 55 sunrise photos taken between June 22 and December 21, measuring the distance of sunrise in each photo from that on the equinox (point D, March 22).

## VIII. COMMENTS AND CONCLUSIONS

When I began this project, my first calculations predicted 8.4 h of daylight for December 21, while the local newspaper reported a time of 8.63 h. This suggested that the simplified model of the Earth-Sun system presented here, and the introduction of the light-dark circle (LDC), could provide a straightforward explanation of the seasonal variations of sunlight. Using rotation matrices rather than vector analysis, diagrams such as Fig. 5 were generated to explain the duration of daylight in the northern hemisphere throughout the year. I have not seen diagrams similar to Fig. 5 in any textbook, so perhaps they introduce a new way to *visualize* seasonal phenomena. Further analysis using rotation matrices allows one to explain the observed path of the Sun from any location on Earth, and the results are in good agreement with sunrise photographs taken throughout the year.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts of interest to disclose.

### APPENDIX A: ROTATION MATRICES

**has known components (**

*r**x*

_{B},

*y*

_{B,}

*z*

_{B}); and system C, with dashed lines, in which

**has**

*r**unknown*components (

*x*

_{C},

*y*

_{C},

*z*

_{C}). In Fig. 11(a), system C is obtained by rotating the X

_{B}and Z

_{B}axes counterclockwise by a (positive) angle

*θ*. Because the Y

_{B}and Y

_{C}axes are parallel, $ y B= y C$. The components of

*z*

_{B}along the X

_{C}and Z

_{C}axes, shown with thick dashed lines, are

*x*

_{C}= −sin

*θ z*

_{B}and

*z*

_{C}= cos

*θ z*

_{B}. The components of

*x*

_{B}along the X

_{C}and Z

_{C}axes, shown with prominent dotted lines, are

*x*

_{C}= cos

*θ x*

_{B}and

*z*

_{C}= sin

*θ x*

_{B}. Summing the components on the X

_{C}and Z

_{C}axes, we find the following:

*r*_{C}and

*r*_{B}are column vectors and

**R**

_{CB}is the 3 × 3 rotation matrix.

*θ*about the X

_{B}axis. The results are given in matrix form by

*θ*about the Z

_{B}axis, shown in Fig. 11(c), are given by

## REFERENCES

*Practical Astronomy with Your Calculator or Spreadsheet*

*Celestial Calculations: A Gentle Introduction to Computational Astronomy*

*Mathematical Methods in the Physical Sciences*

*Excel*” for information about error analysis.