{"id":10243,"date":"2025-12-05T12:25:24","date_gmt":"2025-12-05T17:25:24","guid":{"rendered":"https:\/\/blogs.mathworks.com\/community\/?p=10243"},"modified":"2025-12-06T11:40:05","modified_gmt":"2025-12-06T16:40:05","slug":"obliquity","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/community\/2025\/12\/05\/obliquity\/","title":{"rendered":"The Perils of Obliquity"},"content":{"rendered":"
\n
You might think you don't care about spherical geometry, but really you do. Because spherical geometry is where you live.<\/div>\n
This is a story of the real-world complications of three dimensional projections from one great circle to another. It's also a good example of why I love MATLAB. I can puzzle about an interesting problem, then use my favorite technical computing tool to come to a deeper understanding. I'm going to bring you along on a ride with me as we explore something that I think about every December: the mysteries of the shortest day. Along the way, we'll savor some MATLAB morsels. Lucky you! Buckle up.<\/div>\n
Here's a picture of the sun and the earth at the solstice.<\/div>\n
\n
\n
t = linspace(0,2*pi);<\/span><\/div>\n<\/div>\n
\n
ct = cos(t);<\/span><\/div>\n<\/div>\n
\n
st = sin(t);<\/span><\/div>\n<\/div>\n
\n
<\/div>\n<\/div>\n
\n
clf<\/span><\/div>\n<\/div>\n
\n
heq = line(ct,st,0*t,Color=\"black\"<\/span>);<\/span><\/div>\n<\/div>\n
\n
earth = [0 0];<\/span><\/div>\n<\/div>\n
\n
sun = [-1 0];<\/span><\/div>\n<\/div>\n
\n
line( ...<\/span><\/span><\/div>\n<\/div>\n
\n
[earth(1) sun(1)],[earth(2) sun(2)], ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=0.8*[1 1 1])<\/span><\/div>\n<\/div>\n
\n
he = line(earth(1),earth(2), ...<\/span><\/span><\/div>\n<\/div>\n
\n
Marker=\".\"<\/span>,MarkerSize=24, ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=[0.3 0.6 1]);<\/span><\/div>\n<\/div>\n
\n
hs = line(sun(1),sun(2), ...<\/span><\/span><\/div>\n<\/div>\n
\n
Marker=\".\"<\/span>,MarkerSize=24, ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=[0.95 0.8 0.3]);<\/span><\/div>\n<\/div>\n
\n
axis equal<\/span><\/span><\/div>\n<\/div>\n
\n
axis off<\/span><\/span><\/div>\n<\/div>\n
\n
legend([he,hs,heq],{\"Earth\"<\/span>,\"Sun\"<\/span>,\"Ecliptic\"<\/span>})<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
After one spin on its axis, twenty-four hours later, the sun is here relative to the earth. A day has passed.<\/div>\n
\n
\n
t2 = linspace(pi,1.5*pi,5);<\/span><\/div>\n<\/div>\n
\n
ct2 = cos(t2);<\/span><\/div>\n<\/div>\n
\n
st2 = sin(t2);<\/span><\/div>\n<\/div>\n
\n
sun = [ct2(2) st2(2)];<\/span><\/div>\n<\/div>\n
\n
line([earth(1) sun(1)],[earth(2) sun(2)], ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=0.8*[1 1 1]);<\/span><\/div>\n<\/div>\n
\n
set(hs,XData=sun(1),YData=sun(2))<\/span><\/div>\n<\/div>\n
\n
legend([he,hs,heq],{\"Earth\"<\/span>,\"Sun\"<\/span>,\"Ecliptic\"<\/span>})<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
Okay, I admit I'm making it look like the sun goes around the earth, but relatively speaking this is true enough. The sun icon is indicating the direction of the sun from earth's point of view. This geocentric view makes the diagrams simpler. Also, I'm making the days abnormally long. Pretend there are exactly 16 days in year. So there are only 4 days between the solstice and the equinox.<\/div>\n
The plane in which the earth moves around the sun (represented by the black circle we see here) is called the ecliptic.<\/div>\n
\n
\n
for <\/span>n = 1:length(t2)<\/span><\/div>\n<\/div>\n
\n
line([earth(1) ct2(n)],[earth(2) st2(n)], ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=0.8*[1 1 1])<\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n<\/div>\n
\n
sun = [ct2(end) st2(end)];<\/span><\/div>\n<\/div>\n
\n
set(hs,XData=sun(1),YData=sun(2))<\/span><\/div>\n<\/div>\n
\n
legend([he,hs],{\"Earth\"<\/span>,\"Sun\"<\/span>})<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
The earth rotating on its axis determines how long an hour is. This is the \"clock plane\" (also known as the equatorial plane). The earth's motion around the sun determines how long a year is. This is the \"year plane\" (the ecliptic plane). In the diagram below, the outer circle is the ecliptic plane and we've added a gray inner circle for the equatorial plane.<\/div>\n
As you know, the earth is tilted relative to the ecliptic plane. This angle between the equatorial plane and the ecliptic plane, called obliquity, is approximately 23.4 degrees.<\/div>\n
But first let's consider how things would look if there was no axis tilt. The equatorial plane would be coincident with the ecliptic. There would be no seasons. No solstice. No equinox (or rather, never-ending equinox). A little boring, maybe, but the astronomical calculations would be a breeze. The days marked off on the equatorial plane perfectly match the same days marked off on the ecliptic.<\/div>\n
\n
\n
clf<\/span><\/div>\n<\/div>\n
\n
hecl = line(ct,st,0*t, ...<\/span><\/span><\/div>\n<\/div>\n
\n
Color=[0 0 0]);<\/span><\/div>\n<\/div>\n
\n
for <\/span>n = 1:length(t2)<\/span><\/div>\n<\/div>\n
\n
line(ct2(n)*[0 1],st2(n)*[0 1],[0 0],Color=0.8*[1 1 1])<\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n<\/div>\n
\n
<\/div>\n<\/div>\n
\n
hgt = hgtransform;<\/span><\/div>\n<\/div>\n
\n
heq = patch(0.5*ct,0.5*st,0*t, ...<\/span><\/span><\/div>\n<\/div>\n
\n
FaceColor=0.95*[1 1 1],EdgeColor=[0 0 0], ...<\/span><\/span><\/div>\n<\/div>\n
\n
Parent=hgt);<\/span><\/div>\n<\/div>\n
\n
hpole = line([0 0],[0 0],[0 1],Parent=hgt);<\/span><\/div>\n<\/div>\n
\n
for <\/span>n = 1:length(t2)<\/span><\/div>\n<\/div>\n
\n
line(ct2(n)*[0 0.5],st2(n)*[0 0.5],[0 0],Color=0.8*[1 1 1],Parent=hgt)<\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n<\/div>\n
\n
<\/div>\n<\/div>\n
\n
axis equal<\/span><\/span><\/div>\n<\/div>\n
\n
axis off<\/span><\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
But now let's add the axis tilt to the equatorial plane. I'm going to exaggerate the tilt here to emphasize the discrepancy.<\/div>\n
Oh, and watch how I use HGTRANSFORM<\/a> to rotate a bunch of objects in 3-space in one slick maneuver! HGTRANSFORM is a favorite of mine.<\/div>\n
The key thing to notice here is that the gray \"noon\" lines no longer match up. That means 24 hours (one notch on the inner circle) doesn't match the time from solar noon to solar noon (one notch on the outer circle). The day length is literally stretching and shrinking! You can also see that the lines do match on the solstices and the equinoxes.<\/div>\n
\n
\n
epsilon = 45; % Actual epsilon is only 23.4 deg<\/span><\/span><\/div>\n<\/div>\n
\n
rot = makehgtform(YRotate=deg2rad(epsilon));<\/span><\/div>\n<\/div>\n
\n
set(hgt,Matrix=rot)<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
It's a little hard to visualize what's going on from this top view. Let's look at things from the side.<\/div>\n
\n
\n
view(3)<\/span><\/div>\n<\/div>\n
\n
view([15 20])<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
This is the punchline: we have two competing ways to define a day: the time between successive noons, which varies with the season, and exactly 24 hours, which is constant throughout the year. Things look good at the solstice and the equinox, but significantly off in between.<\/div>\n
This is a real effect. Solar days vary in length. Obliquity contributes to the so-called Equation of Time<\/a> which spells out when the sun is running ahead of or behind what the clock says.<\/div>\n
The code below will generate an animation. I'll paste in the corresponding GIF-file.<\/div>\n
\n
\n
for <\/span>epsilon = 45*sin(linspace(0,pi))<\/span><\/div>\n<\/div>\n
\n
rot = makehgtform(YRotate = deg2rad(epsilon));<\/span><\/div>\n<\/div>\n
\n
set(hgt,Matrix=rot)<\/span><\/div>\n<\/div>\n
\n
axis([-1 1 -1 1 -1 1])<\/span><\/div>\n<\/div>\n
\n
view(2)<\/span><\/div>\n<\/div>\n
\n
drawnow<\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
\"\"<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
How bad is the discrepancy? Up to 10 minutes at its maximum. Let's try a couple of alternative tilt scenarios. If the earth was tilted ten degrees more, the solar noon would be up to 20 minutes out of sync with the clock.<\/div>\n
\n
\n
clf<\/span><\/div>\n<\/div>\n
\n
<\/div>\n<\/div>\n
\n
t = linspace(pi,3\/2*pi);<\/span><\/div>\n<\/div>\n
\n
discrepancy = zeros(size(t));<\/span><\/div>\n<\/div>\n
\n
for <\/span>epsilon = [13.4 23.4 33.4]<\/span><\/div>\n<\/div>\n
\n
m = makehgtform(YRotate = deg2rad(epsilon));<\/span><\/div>\n<\/div>\n
\n
<\/div>\n<\/div>\n
\n
for <\/span>n = 1:length(t)<\/span><\/div>\n<\/div>\n
\n
p1 = [cos(t(n)); sin(t(n)); 0; 1];<\/span><\/div>\n<\/div>\n
\n
p2 = m*p1;<\/span><\/div>\n<\/div>\n
\n
% Dropping the z term effectively projects the vector onto the x-y plane<\/span><\/span><\/div>\n<\/div>\n
\n
p2(3) = 0;<\/span><\/div>\n<\/div>\n
\n
% Now find the angle between the two vectors<\/span><\/span><\/div>\n<\/div>\n
\n
v1 = p1(1:3);<\/span><\/div>\n<\/div>\n
\n
v2 = p2(1:3);<\/span><\/div>\n<\/div>\n
\n
discrepancy(n) = atan2(norm(cross(v1,v2)), dot(v1,v2));<\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n<\/div>\n
\n
plot(t,discrepancy\/(2*pi)*24*60)<\/span><\/div>\n<\/div>\n
\n
hold on<\/span><\/span><\/div>\n<\/div>\n
\n
end<\/span><\/span><\/div>\n<\/div>\n
\n
set(gca,XLim=pi*[1 1.5], ...<\/span><\/span><\/div>\n<\/div>\n
\n
XTick=pi*[1 1.25 1.5], ...<\/span><\/span><\/div>\n<\/div>\n
\n
XTickLabels={\"Solstice\"<\/span>,\"Cross-Quarter Day\"<\/span>,\"Equinox\"<\/span>})<\/span><\/div>\n<\/div>\n
\n
ylabel('Discrepancy (minutes)'<\/span>);<\/span><\/div>\n<\/div>\n
\n
title('Time Discrepancy'<\/span>);<\/span><\/div>\n<\/div>\n
\n
grid on<\/span><\/span><\/div>\n<\/div>\n
\n
legend({'\\epsilon = 13.4 deg'<\/span>,'\\epsilon = 23.4 deg'<\/span>,'\\epsilon = 33.4 deg'<\/span>})<\/span><\/div>\n
\n
\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n
<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
What's going on here? One hour is NOT 1\/24 of the time between successive noons. If it was, we'd go insane. To keep from having an hour that is constantly changing, humans created an invariant average \"standard hour\". The obliquity effect stretches and shrinks the solar day as measured by this standard hour. This moves the local solar noon around, depending on the time of year. And one effect of this is that the earliest sunset arrives sooner than the shortest day. And in fact, the earliest sunset at my latitude is December 9th.<\/div>\n
Happy Crepusculus<\/a>!<\/div>\n<\/div>\n