Greetings from Natick, Massachusetts, where the outdoor temperature is currently a balmy 15 degrees Fahrenheit and the sun sets at 4:26 in the afternoon. The days are getting colder but at least they are getting longer now. Tomorrow the sun will set a minute later. Tomorrow’s sunrise, on the other hand, will be at more or less the same time as today’s. In fact, at this latitude, today is the latest sunrise of the year: 7:15 AM.

These short days always make me think about the geometry of earth and sun that gives rise to seasons. And this relates to MATLAB because it gave me the opportunity to answer a question and write a little code. The question is this: for any given day, let’s say today, where are the people who will see the sun rise at exactly the same moment as me? Answering this question gave me the chance to learn a little Mapping Toolbox code.

First some basic geometry. The earth is always half dark and half illuminated by the sun. So the sunset/sunrise line (sometimes called the gray line) that separates day from night is a great circle, which is to say a circle that divides the earth evenly into two halves. This gray line is centered on the point on the earth that is precisely under the sun right now. Since we are near the winter solstice (for the northern hemisphere), the sun can be found hovering close to the Tropic of Capricorn at around 22.8 degrees south latitude. So now we have a goal: draw a great circle with its center positioned at -22.8 deg. latitude and a longitude such that the sunrise line passes through Natick.

Now let’s learn some basic Mapping Toolbox skills. This was fun for me, because two hours ago I didn’t know anything about the Mapping Toolbox. After a quick search for great circles, I started with an example from a page on circles in the documentation. From there I quickly got a map like this.

clf
axesm('MapProjection','ortho')
gridm on
framem on
setm(gca, ...
    'Origin', [30 -50 20], ...
    'MLineLimit', [75 -75], ...
    'MLineException',[0 90 180 270])

landareas = shaperead('landareas.shp','UseGeoCoords',true);
geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);

Let’s put Natick on the map. That’s me under the red dot.

natick = [42.3 -71.4];
plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')

Here’s something really cool about the Mapping Toolbox: I can do a completely different map projection by changing one line in the above code. Here is the Miller Cylindrical projection. Fun!

clf
axesm('MapProjection','miller')
gridm on
framem on
setm(gca, ...
    'Origin', [30 -50 20], ...
    'MLineLimit', [75 -75], ...
    'MLineException',[0 90 180 270])
geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);
plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')

I could also do a Murdoch III Minimum Error Conic Projection, but that would just be showing off.

Now down to business. Sunrise in Natick is at 7:15 Eastern Standard Time. At that moment, the sun is above the south Atlantic at 22.8 deg. south latitude and 4 deg. west longitude. We’ll use the function SCIRCLE1 to generate the great circle.

clf
axesm('MapProjection','ortho')
gridm on
framem on
setm(gca, ...
    'Origin', [30 -50 20], ...
    'MLineLimit', [75 -75], ...
    'MLineException',[0 90 180 270])
sunLocation = [-23 -4];
grayLine = scircle1(sunLocation(1), sunLocation(2), 90);

geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);

plotm(sunLocation(1), sunLocation(2),'ko','MarkerFaceColor','y')
plotm(grayLine(:,1), grayLine(:,2),'red')
plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')

Now we’re getting somewhere! Let’s switch back to the Miller projection and zoom in to the east coast of the U.S.

clf
axesm('MapProjection','Miller', ...
    'MapLatLimit',[-11 60], ...
    'MapLonLimit',[-50 40])
gridm on
framem on
setm(gca, ...
    'Origin', [11 -59 8], ...
    'MLineLimit', [75 -75], ...
    'MLineException',[0 90 180 270])
sunLocation = [-23 -4];
grayLine = scircle1(sunLocation(1), sunLocation(2), 90);

geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);

plotm(sunLocation(1), sunLocation(2),'ko','MarkerFaceColor','y')
plotm(grayLine(:,1), grayLine(:,2),'red')
plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')

And there’s the answer we’ve been after. Tomorrow, in the unlikely event I get up early enough to see the sun rise, I’ll be sharing that experience and that instant with people in Arroyos De Mantua, Cuba and Ísafjörður, Iceland. There’s something comforting about that.

Published with MATLAB® R2013b