I enjoy astronomy. In this post I’m going to show off a cool file from the File Exchange and demonstrate MATLAB’s new-ish datetime functionality (including dateshift and duration). But it’s also a fun opportunity to tell a little astronomical story. A story that starts with a simple question: Who gets to decide when it’s noon? The sun, or your wristwatch?
Noon is defined to be the moment that divides morning from afternoon. Astronomically speaking, this is when the sun stops getting higher and starts getting lower. But practically speaking, it’s when Mickey’s hands point to twelve.
Who is right?
Mickey is right. How can that be? Why should the watch on your wrist have more authority than the star that we orbit?
Consider that if your timepiece is a sundial, there’s no conflict. By definition, it’s noon when the sun reaches its highest point in the sky, or when the gnomon’s shadow is the shortest (the gnomon is the pointy thing that casts the shadow). Sundials are simple and sturdy, but they have some issues. They don’t work very well at night, but more to the point, they imply an infinite number of time zones, one for every possible line of longitude.
Even this inconvenience didn’t really matter until the arrival of the railroad and the telegraph. Only then, and for the first time in history, did people travel so quickly that anyone actually cared that noon in New York wasn’t quite the same as noon in Washington DC. In particular, constructing railroad timetables was a nightmare.
To avoid the infinite time zones problem, we bin meridians into 24 groups, one for each hour of the day. One special meridian in each time zone uses the sun to determine noon, and the rest of the zone falls into line. Problem solved, right? Not quite.
Here’s where we need to use some MATLAB. I’ll be using Darin Koblick’s excellent Vectorized Solar Azimuth and Elevation Estimation. I added the file using the Add-On Explorer. Here’s what it looks like in the Add-On Manager.
This file can tell me where the sun is right now. It returns the answer as azimuth (north = 0 deg, east = 90 deg) and elevation (horizon = 0 deg, straight up = 90 deg).
Here are the coordinates where I sit.
lat = 42.299899; % deg lon = -71.350729; % deg alt = 0.305; % km tz = hours(timezone(lon)); % time zone adjustment from universal time
At this instant, where is the sun?
dtNow = datetime('now'); dtNow.Format = 'yyyy/MM/dd HH:mm:ss'; [az,el] = SolarAzEl(dtNow+tz,lat,lon,alt); fprintf('Solar Azimuth: %5.2f deg\nSolar Elevation: %5.2f deg\n',az,el)
Solar Azimuth: 221.82 deg Solar Elevation: 17.26 deg
It’s a little west of south (mid-afternoon) and somewhat low in the sky, as befits late January.
Here is the elevation of the sun in the sky for every minute between 11 AM and 1 PM local time. By taking the maximum elevation and finding the corresponding time, we’ve determined local noon, which is to say, the astronomical moment at which the sun is highest.
% Make a list of all the minutes from 11 AM to 1 PM timeList = datetime(year(dtNow),month(dtNow),day(dtNow),11,0:120,0)'; % Find the corresponding solar elevations [az,el] = SolarAzEl(timeList+tz,lat,lon,alt); sun = [0.96 0.74 0.24]; sky = [0.85 1 1]; plot(timeList,el, ... 'LineWidth',1,'Color',sun) set(gca,'Color',sky) [~,ix] = max(el); noon = timeList(ix); noonEl = el(ix); line(noon,noonEl, ... 'Color',sun, ... 'Marker','.','MarkerSize',40) grid on xlabel('Time') ylabel('Elevation Above Horizon (deg)') noon.Format = 'HH:mm'; title(sprintf('Local Noon Occurs at %s',noon))
It’s close, but “noon” isn’t quite noon. Part of that effect is because of where I sit in the time zone. But that’s not all there is to it.
Let’s look at what happens over a period of few days. These are the solar elevations for the first two weeks of the year. You can see noon is shifting slightly later every day, and the sun is climbing slightly higher every day.
noons = datetime(year(dtNow),1,2:14); for n = 1:length(noons) timeList = datetime(year(noons(n)),month(noons(n)),day(noons(n)),11,30:90,0)'; [~,el] = SolarAzEl(timeList+tz,lat,lon,alt); [~,ix] = max(el); plot(timeList-dateshift(timeList,'start','day'),el,'Color',sun) hold all plot(timeList(ix) - dateshift(timeList(ix),'start','day'),el(ix), ... 'Marker','.','MarkerSize',32, ... 'Color',sun) end hold off grid on xlabel('Time') ylabel('Elevation Above Horizon (deg)') set(gca,'Color',sky)
It looks a little bumpy because the times are discretized to the nearest minute.
Now let’s plot where noon occurs for the entire year.
noons = datetime(year(dtNow),1,1:7:365); els = zeros(size(noons)); for n = 1:length(noons) timeList = datetime(year(noons(n)),month(noons(n)),day(noons(n)),11,0:120,0)'; [~,el] = SolarAzEl(timeList+tz,lat,lon,alt); [~,ix] = max(el); noons(n) = timeList(ix); els(n) = el(ix); end plot(noons - dateshift(noons,'start','day'),els, ... 'Marker','.','MarkerSize',32, ... 'Color',sun) set(gca,'Color',sky) grid on xlim(duration(hours([11 12.5]))) xlabel('Time') ylabel('Elevation (degrees)') title('Elevation of the Sun at Local Noon')
Yow! Where did that figure 8 shape come from?
mean(noons - dateshift(noons,'start','day'))
ans = duration 11:45:29
On average, local noon occurred at 11:45. But sometimes it was earlier and sometimes it was later by as much as 15 minutes. There’s clearly more going on here than just time zone shifts.
There are two subtle effects here. One is due to the fact that the earth’s orbit isn’t quite circular. It’s slightly elliptical, so the earth actually moves faster around the sun at certain times of the year. The other effect is due to the fact that the earth is tilted on its axis by 23.4 deg. These two factors make a the length of a day, as measured from local noon to local noon, vary by as much as a half hour. Crazy, right?
At some point hours and minutes were redefined as fractions of a single standard “mean solar day” rather than a locally measured (but inconsistent) “sundial day”. We’ve ironed out all the wrinkles in time. Put another way, first the sun defined minutes, then minutes defined the sun. The timekeeping tail wags the dog.
This is pretty basic stuff, as far as astronomy goes, but it was so much fun to find a nice file on the File Exchange and then be able to bang out this code so quickly. It’s amazing how much good stuff is waiting for you there. Go look!
4 CommentsOldest to Newest
Great post! I see on the FEX that this algorithm is accurate to +/- 1 degree. Those interested in calculating solar zenith angles more accurately can look here: http://rredc.nrel.gov/solar/codesandalgorithms/spa. The authors demonstrate an accuracy of +/- 0.0003 degrees. They provide C-code that can be converted to Matlab or a MEX file.
Note that in their PDF, the value for H in Table A5.1 of Section A5 has a typo. It should be 11.105902 degrees rather than 11.105900. This is corrected in their code and is necessary to obtain their answers in the examples.
Interestingly, I also found that the algorithm in their paper for calculating Julian dates disagrees with Matlab’s juliandate() function for dates before October 15, 1582. The discrepancy has to do with the 10 day error in the Julian calendar in 1582 when the Gregorian correction was made, I believe. See this PNG on Wikipedia.
Thanks for the comments Eric! In astronomy there are always opportunities for more precision. Like your Gregorian calendar story, I’ve always been amused about how astronomers and historians differ about dates before 1 CE. For historians, there is no year zero. 1 BCE gives way to 1 CE. Astronomers, on the other hand, insist on using math with a zero between -1 and 1! So you have to apply a correction factor to make the two square with one another.
Great article. I was reminded of this other article that illustrates deviations in solar noon across the world:
I love that time-wrongness map. China really jumps out for being one gigantic time zone. And you can see the world is mostly skewed late (red) rather than early (green). Thanks for sharing!