Daytime length variations with latitude and season (1/3)

[This article was first published on R on BitFoam, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

It is easy to explain the changing seasons and the variation in the length of day and night time at different latitudes by means of the axial tilt, also known as obliquity (Wikipedia (n.d.a)).

Earth currently has a mean tilt of 23.44° (Wikipedia (n.d.a)) between:

  1. its rotational axis, the imaginary line that passes through both poles, and on which Earth spins producing the 24-hour day,
  2. and its orbital axis, the line perpendicular to the imaginary plane (ecliptic) across which Earth moves as it revolves around the Sun, producing the 365-day year.

The figure 1 shows the two axes and the axial tilt between them.

Axial Tilt, By I, Dennis Nilsson, CC BY 3.0, https://commons.wikimedia.org/w/index.php?curid=3262268

Figure 1: Axial Tilt, By I, Dennis Nilsson, CC BY 3.0, https://commons.wikimedia.org/w/index.php?curid=3262268

Over the course of a revolution around the Sun, the orientation of the rotational axis remains the same respect to the background of distant stars. As a result, the orientation of the rotaxional axis relative to the Sun varies along the Earth’s orbit.

The figure 2 shows the orientation of the rotational axis in two opposite points of the Earth’s orbit, on June solstice and December solstice, when the tilt of the rotational axis is at its maximum, toward the Sun or away from the Sun.

Axis Orientation in opposite points of the Earth's orbit, By cmglee, NASA - Own work using: http://visibleearth.nasa.gov/view.php?id=73580, Public Domain, https://commons.wikimedia.org/w/index.php?curid=41309095

Figure 2: Axis Orientation in opposite points of the Earth’s orbit, By cmglee, NASA – Own work using: http://visibleearth.nasa.gov/view.php?id=73580, Public Domain, https://commons.wikimedia.org/w/index.php?curid=41309095

This post is the first of a series that explores how to calculate and visualize the daytime length by latitude and day in the year using R.

Simple model

A simple model has been built upon the concepts described above. The most relevant approximations in the model are listed below together with a short note qualifying to what extent deviates from the reality.

  1. Earth is modelled like a sphere. In fact, Earth is only approximately spherical, it deviates mainly due to the centrifugal force caused by its rotation (Wikipedia, n.d.c) and, to a lesser extent, to its non-uniform density and geological phenomena.
  2. Earth spins on its own axis (the rotational axis) once every 24 hours. Earth’s rotational period is not constant mainly due to the tidal forces of the Moon (Wikipedia, n.d.b).
  3. Earth revolves about the Sun in a circular orbit once every 365 days. Celestial orbits are not perfectly circular. While Earth’s orbit is better approximated by an ellipse, the influence of the other solar system bodies continously alters its orbit.
  4. Earth’s motion is decomposed in an intertwined succession of discrete rotations and translations (revolution): rotates, revolves, rotates, revolves … Earth’s motion is continous, it spins as it revolves.
  5. The axial tilt is 23.44° and remains the same over an orbital period. The axial tilt oscillates between 22.1° and 24.5° on a 41,000-year cycle (Wikipedia (n.d.a)).
  6. Parallel sunlight. The Sun and the Earth are so far apart that the sunlight appears to be parallel on Earth. Check this animation to get a feeling of how good is this approximation.
  7. At any time, half of the Earth is illuminated by the Sun and the other half is at night. This removes the necessity of modelling the Sun at all. This is a result of simplifications 1) and 6).
  8. The effect of the atmospheric refraction is negligible in the day- and nightime. The refraction may make the sun visible for several minutes before it raises and after it sets (Konstantin Bikos, n.d.).

Conventions

The figure 3 depicts the conventions in relation with the system coordinates followed in the series.

  • x-axis and y-axis make up the horizontal plane
  • z-axis is always the vertical axis and is positive upwards.
  • The polar coordinate \(\theta\) is measured from the z-axis (positive pointing upwards).
  • Rotation follows the right-hand rule, when the thumb is pointing in the positive direction along the axis of rotation, then the curled fingers indicate the positive rotation direction.
Cartesian and polar coordinates, By Andeggs - Own work, Public Domain, https://commons.wikimedia.org/w/index.php?curid=7478049

Figure 3: Cartesian and polar coordinates, By Andeggs – Own work, Public Domain, https://commons.wikimedia.org/w/index.php?curid=7478049

Grid generation

The first step is the generation of points to model Earth, these points are then used to determine the daytime length at different latitudes at different times in the year.

The grid is determined by the following parameters:

Latitude points

Set of angles raging from -90° (South Pole) to 90° (North Pole), 0° at the Equator. The number of latitude points is controlled by the parameter lat_step. The parameter is set to 10, resulting in 19 latitudes, to keep a manageable number of points and a resonable time for rendering the plots. The value of the step can be reduced to increase the number of points and the resolution.

1lat_step <- 10
2lat_pts <- seq(from = - 90, to = 90, by = lat_step)
3num_lat_pts <- length(lat_pts)

Calculating latitude points

The polar coordinate \(\theta\) can be determined from the latitude in degrees using the following formula: \(\theta (rad) = {\pi \over 2} \times { (1 - {latidude (degrees) \over 90})}\)

Once we have the angle \(\theta\), the cartesian coordinates may be retrieved from the spherical coordinates by the following equations:

  • \(x = r \times cos(\varphi) \times sin(\theta)\)
  • \(y = r \times sin(\varphi) \times sin(\theta)\)
  • \(z = r \times cos(\theta)\)

Thus, for \(r = 1\)

1coor_lat <- matrix(1, num_lat_pts, 3)
2x_lat <- sin( pi / 2 * (1 - lat_pts / 90))
3y_lat <- sin( pi / 2 * (1 - lat_pts / 90))
4z_lat <- cos( pi / 2 * (1 - lat_pts / 90)) 
5coor_lat[,1] <- x_lat
6coor_lat[,2] <- y_lat
7coor_lat[,3] <- z_lat 
8colnames(coor_lat) <- c('x', 'y', 'z')

… and \(\varphi = 0\) the figure 4 plots in 3D the resulting latitude points using the library plotly.

  • \(x = r \times sin(\theta)\)
  • \(y = 0\)
  • \(z = r \times cos(\theta)\)
1scene = list(camera = list(eye = list(x = 0, y = 2, z = 0)))
2fig_lat <- plot_ly() %>% config(displayModeBar = FALSE, scrollZoom = FALSE)
3fig_lat <- fig_lat %>% add_trace(data = coor_lat, x = coor_lat[,'x'], y = 0, z = coor_lat[,'z'], type ='catter3d', mode='markers', size=1) %>% layout(scene = scene)
4fig_lat

Figure 4: Latitude points

Points along the parallel (longitude)

Set of angles into which a parallel (line of constant latitude) circling the Earth is divided. To simplify the calculation of the daytime length later in this series it ranges from 1 to 48 (multiple of 24).

The number of points is controlled by the parameter lon_step. The parameter is set to 1, resulting in 48 points for each parallel.

1lon_step <- 1
2lon_pts <- seq(from = 1, to = 48, by = lon_step)
3num_lon_pts <- length(lon_pts)

Calculating points on the Earth’s surface

For each latitude point a parallel line consisting of num_lon_pts is created. Therefore Earth´s surface is represented by num_lat_pts * num_lon_pts points of its spherical surface.

In order to do this, each latitude point is repeated num_lon_pts times. The \(cos(\varphi)\) and \(sin(\varphi)\) factors required to build a cicle are generated from the set of points along the parallel and applied to x and y coordinates for each latitude according to:

  • \(x = r \times cos(\varphi) \times sin(\theta)\)
  • \(y = r \times sin(\varphi) \times sin(\theta)\)

This code shows a very efficient way of doing it by using the function rep.

1coor_lat_lon <- matrix( rep( coor_lat, each = num_lon_pts), (num_lat_pts * num_lon_pts), 3) 
2colnames(coor_lat_lon) <- c('x', 'y', 'z') 
3coor_lat_lon[,'x'] <- coor_lat_lon[,'x'] * rep( cos( lon_pts * 2 * pi / num_lon_pts), num_lat_pts) 
4coor_lat_lon[,'y'] <- coor_lat_lon[,'y'] * rep( sin( lon_pts * 2 * pi / num_lon_pts), num_lat_pts) 

The figure 5 plots the resulting points.

1scene = list(camera = list(eye = list(x = 0, y = 2, z = 0))) 
2fig_lat_lon <- plot_ly() %>% config(displayModeBar = FALSE, scrollZoom = FALSE) 
3fig_lat_lon <- fig_lat_lon %>% add_trace(data = coor_lat_lon, x = coor_lat_lon[,'x'], y = coor_lat_lon[,'y'], z = coor_lat_lon[,'z'], type ='scatter3d', mode='markers', size=1) %>% layout(scene = scene) 
4fig_lat_lon 
To leave a comment for the author, please follow the link and comment on their blog: R on BitFoam.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)