The Longest Day of the Year

The following description will walk you through the code on my GitHub page for calculating the differences in the length of a day on Earth as a function of distance from the Sun!

The code presented and explained herein describes the motion of the earth as it orbits the Sun, and as it rotates on its own axis. The calculations utilize Kepler’s laws of planetary motion and the principle of conservation of angular momentum to show that the longest day of the year occurs at the beginning of January, approximately two weeks after the winter solstice.

Contrary to popular belief, the longest day of the year is not the summer solstice in June! The summer solstice is the day where we (in the northern hemisphere, at least) experience the greatest amount of daylight. This increased daylight is due to the tilt of the Earth on its axis, which is what is responsible for the changing of seasons through the year. The calculations here pertain, instead, to the day as defined by the rotation of the Earth on its own axis, where one full rotation is equivalent to a day. Therefore, the hours of sunlight are a completely separate, and in this case, unimportant, phenomenon.

According to Kepler’s first law of motion, a planet is defined as a body that orbits its star in an elliptical pattern, with the Sun at one focal point. Below is a representation of what this looks like, with the eccentricity (or, “the degree of ellipticity”) greatly exaggerated for improved visualization:

Exaggerated depiction of the Earth’s orbit about the Sun — the orbit is an ellipse with the Sun at one focal point.

The eccentricity of a shape indicates how “oval” the shape is. A circle has an eccentricity of 0, while a straight line has an eccentricity of 1. All other numbers between 0 and 1 correspond to an ellipse with varying degrees of “oval-ness”. The eccentricity of the Earth is quite small, roughly 0.0167, and therefore, the orbit of the Earth around the Sun may appear circular. However, as we have learned just now, this is not the case.

In the above schematic, we have the Sun at one focal point (as defined by Kepler’s first law of motion), with the red circle representing Earth at perihelion (the point closest to the Sun), and the green circle representing aphelion (the point farthest away from the Sun).

Again, the two types of motion we must concern ourselves with are orbiting and rotation. Orbiting refers to the Earth’s travel around the Sun, while rotation refers to the Earth spinning on its own axis. Because of the elliptical nature of the Earth’s motion about the Sun, the days are not all the same length. This may seem strange, because we on Earth define our days to be exactly 24 hours, every day. So, what causes the days to be different in length, and how different are the lengths of the days?

To begin, we must first remember that the law of conservation of angular momentum states that angular momentum of a closed system will remain constant. Angular momentum can be described by the equation:

Conservation of Angular Momentum

Where m is the mass of the body (in this case, the mass of the Earth),

v is the velocity of the Earth around the Sun (orbital velocity), and

r is the distance between the Earth and the Sun, which we know, due to the elliptical orbit of the Earth around the Sun, is not the same every day.

Because the mass of the Earth is constant, we can ignore it when using the law of conservation of angular momentum. The orbital velocity and distance are not constant, and therefore, we focus our attention on these two variables.

When the Earth is closest to the Sun, at perihelion, it must then have a larger orbital velocity compared to when it is at aphelion, to maintain conservation of angular momentum. This means that the velocity and distance at perihelion multiplied together is equal to the velocity and distance at aphelion multiplied together:

As the Earth gets closer to the Sun, it moves slightly faster (increased orbital velocity), and when it is farther away, it moves slower.

Looking back at the schematic above, you may notice that at perihelion, the angle for one day (purple) is larger than it is at aphelion (orange). According to Kepler’s second law of motion, a planet will sweep out equal areas in equal time. What this means with regards to our schematic is that if both the purple sector and the orange sector represent the same amount of time, then the areas of these sectors are equal. The only thing that is obviously NOT equal are the angles through which they move, due to the difference in distance from the Sun. This is the angle through which the Earth rotates on its own axis. Notice that it is larger at perihelion than at aphelion.

The different values for Earth’s aphelion distance, perihelion distance, associated orbital velocities and more are included in the class called EarthParameters in the code. Many of these values can be found on the NASA website, or elsewhere on the internet. Not all of them are used in this codebase, but having them in this class is valuable for another calculation that someone may want to work on.

The value for angular momentum, which is constant, is defined in the code using values from the EarthParamters class as:

        constant_Vr = EarthParameters.avg_v * EarthParameters.AU

Here, we have used the average orbital velocity of the Earth around the Sun, and the average distance, defined in astronomy as an Astronomical Unit (AU). This value must stay constant from day to day, lest we have a violation of conservation of angular momentum.

We also need to know about the aphelion day of last year, 2022, and the perihelion day for next year, 2024, in order to calculate all the values we need for the Earth’s motion around the Sun in 2023. This information may also be found online, as the purpose of this exercise is not to determine when aphelion or perihelion will occur, but rather to show mathematically that not all days of the year are the same length.

We can then calculate the following:

distance_p_to_a= (EarthParameters.aph - EarthParameters.peri)/183 
#This is the approximate change in distance per day between 2023's perihelion (January 4th) and aphelion (July 6th), which is 183 days.

distance_a_to_p= (EarthParameters.aph - EarthParameters.peri)/180 
#This is the approximate change in distance per day between 2023's aphelion (July 6th) and 2024's perihelion (January 2nd), which is 180 days.

previous_aph_to_current_peri_distance_per_day=(EarthParameters.aph - EarthParameters.peri)/184 
#This is the approximate change in distance per day between 2022's aphelion (July 4th) and 2023's perihelion (January 4th), which is 184 days.

previous_aph_to_current_peri_radians_per_day=pi/184 #the earth will cover pi radians around the ellipse from aphelion of 2022 to perihelion of 2023; this metric gives the average number of radians it travels per day.

We now want to start creating a dataframe (df) and filling in the values. We start by creating three columns, one for the date, one for the Earth’s orbital velocity, and another for the distance between the Earth and the Sun.

df = pd.DataFrame(columns = ["Date", "Earth_Orbital_Velocity", "Distance"])

We first fill in the values for each day of the year in the Date column of the df:

df['Date']=pd.date_range(start='1/1/2023', end='12/31/2023', freq='D') #creating df column that includes all 365 days of the year
df['Date']=pd.to_datetime(df['Date'])
df['Date']=df['Date'].dt.strftime('%Y-%m-%d')

Next, we must input a few values by hand, since we are dealing with a different aphelion day from last year. The indices of 3 and 186 correspond to the index value for the perihelion row and aphelion row, respectively. Then, we fill in the approximate value for the distance of the first three days of the year before perihelion (in 2023, this will be January 1st, 2nd and 3rd).

df['Distance'][3] = EarthParameters.peri #filling in perihelion distance
df['Distance'][186] = EarthParameters.aph #filling in aphelion distance
#here, we fill in the approximate distance from the Sun in meters for the first three days of the year, before perihelion.
df['Distance'][2] = (EarthParameters.peri+previous_aph_to_current_peri_distance_per_day)
df['Distance'][1] = (EarthParameters.peri+2*previous_aph_to_current_peri_distance_per_day)
df['Distance'][0] = (EarthParameters.peri+3*previous_aph_to_current_peri_distance_per_day

Next, we can fill in the distance values for all days between perihelion and aphelion in 2023, and from aphelion to the end of the year (not quite the next perihelion), which can be done with the following code:

#here, we fill in approximate values of how the distances changes per day between the Earth and the Sun.
for i in range(4, 186):
    df['Distance'][i]=df['Distance'][i-1] + distance_p_to_a 
for i in range(187, 365):
    df['Distance'][i]=df['Distance'][i-1] - distance_a_to_p

Going back to the conservation of angular momentum, we know that if we have the distance between the Sun and the Earth, we can calculate the orbital velocity on that date by using the earlier metric we calculated called constant_Vr! We can also calculate the angular velocity, which is given as:

If you do not remember this equation, don’t worry! It should be in any introductory-level mechanics textbook. We can also calculate the radians of rotation per day by taking the angular velocity (which is in radians/second) and multiplying it by how we define a day — 24 hours or 86,400 seconds.

So:

df['Earth_Orbital_Velocity']=constant_Vr/df['Distance'] #this is earth's orbital velocity. Note that this is different than rotation velocity (spinning on its own axis)
df['Angular_Velocity']= df['Earth_Orbital_Velocity']/df['Distance'] #this is the angular velocity in radians/second
df['Radians_of_Rotation']=df['Angular_Velocity']*86400 #radians the Earth rotates through each day based on our definition of a day being 86400 seconds.

Now comes the challenging part — calculating the distance through which the Earth moves on a day-to-day basis. If this were a circular orbit, it would be easy — just divide 2πr, which is the circumference of a circle, by 365! However, the Universe is not always that kind, and her laws state that this orbit is in an ellipse, for which the sector length is much more difficult to calculate.

I have used an online resource for this, written by someone named John Cook (thanks, John!) where he uses python code to calculate the sector length of an ellipse.

I have used his calculation to create the function:

def arc_length_in_meters(T0, T1, a, b):
        f = lambda t: sqrt((a**2*sin(t)**2 + b**2*cos(t)**2))
        return quad(f, T0, T1)

Here, T0 is the starting angle, T1 the ending angle, a is the semi-major axis and b the semi-minor axis (both defined in the EarthParameters class). I have drawn my beginning schematic to show aphelion at 0 degrees, perihelion at 180 degrees (or π radians), and the “average” day, halfway between the two, at 90 degrees or π/2 radians.

We can then calculate the arc length per day, and add this information to its own column in our df:

arc_lengths=[]
for i in df.iloc:
    start=0 #we can start at 0 because the only difference in angle we care about is the radians of rotation. This makes the integral (and the code) simpler!
    end=i.Radians_of_Rotation
    X= arc_length_in_meters(start, end, EarthParameters.semi_maj, EarthParameters.semi_min)
    X=X[0]
    arc_lengths.append(X)
df['Arc_Length']=arc_lengths

Once we know how many meters the Earth travels each day around the Sun (remember, it won’t be the same every day!), we can then calculate how long it will take to do so based on the Earth’s orbital velocity. I have created different columns in my df to show you how much (or, more importantly, how little) these changes are:

df['Rotation_Time']=df['Arc_Length']/df['Earth_Orbital_Velocity']
df['Difference_from_defined_Day_in_Minutes']=((df['Rotation_Time']-86400)/60)
df['Difference_per_hour_in_minutes']=df['Difference_from_defined_Day_in_Minutes']/24
df['Difference_per_minute']=df['Difference_per_hour_in_minutes']/60
df['Difference_per_second']=df['Difference_per_minute']/60

If we look at the ‘Difference_in_defined_Day_in_Minutes’, we see that some of these days are about 24 whole minutes longer than the “average” day, or what we define as the average day — 24 hours. But if we break it down further, we see that each second of the day is only elongated by 0.0002777 seconds, or 0.2 milliseconds. This is not nearly enough to affect anything significantly, timewise. For example, Olympic running times won’t be affected based on which day you run your race. These times are so small that they won’t be noticeable at all.

One thing I want to note is that there is not a lot of data on exactly how much the days of the earth change on a day-to-day basis, so I realize that my calculations could be off. I felt that these numbers were a little too large to be correct, but if that is the case, that means the effects are even more miniscule than I have shown here. More importantly, I have demonstrated that there are certain days of the year that are longer or shorter than others, which was the point of this exercise, and have created a dataframe that may be of use to someone out there. I hope that this code and explanation were understandable, and thanks for reading! 😊