Is it dark outside? Draw a sun map!

Our closest star, the sun, is quite fidgety. The times it rises and sets depend on where you are, and whether it is winter or not.

We would like to be able to deduce if the sun is shining outside without having to leave the comforts of our basements - which is why we need an up-to-date sun map (a.k.a. daylight map). You are the one to write a program that generates just that!

Rules: Your program should output an image (in a known format) or an ASCII art representation of our planet, showing (an approximation of) which parts are currently lit by the sun. Your program must be original and self-contained: you are not allowed to copy, use, include or call any code except your programming language's standard libraries.

If you still don't have a clue what I'm talking about, here's an example from Wikipedia:

This is a popularity contest. You should note in your answer which of the following you're trying to achieve (multiple choices possible):

• Correctness. Note that the rules say 'an approximation of' - the better your approximation, the more points in this category. You can check your implementation against Wolfram Alpha's, Time and Date's or die.net's.

• Functionality. For example, what about interactivity? Marking specific locations? Mapping other planets?

• Aesthetics. Drawing continents? Bonus points. Textured continents? Bonus points. On a 3D earth? With clouds? Stars? Correct stars? Massive bonus points. And so on.

• Using uncommon, old, or just plain wrong technology. Sure, you could whip this up in Mathematica, but have you considered using m4? SQL? Forth? x86 assembly?

• Fun. Want to use a Dymaxion-projection map? Go ahead!

• Short code. This is the Code Golf SE, after all.

Have fun!

• @PeterTaylor Probably! I think you don't even have to do lighting: you could draw a 3D earth and rotate it so that the daylight part (and nothing else) faces the viewer. It wouldn't show the nighttime part of the planet, but that's not required. Mar 16, 2014 at 18:17
• An alternative solution is simply to install windows. (In the basement, I mean.) Mar 16, 2014 at 21:55
• @qwr I'd take x86 assembly over m4 any day of the week for this task... Mar 17, 2014 at 0:52
• Also: "the sun is quite fidgety because it rises and sets at different times." Definitely the sun's fault :P
– qwr
Mar 17, 2014 at 1:08
• @qwr Ever thought about how easy time/date calculations would be if we had 10-hour days, 10-day weeks, 10-week months and 10-month years, and sun between t=0 and t=5? But no, the sun has to go and show its ugly face at different places at different times and take too long to go around the earth. No complaints division either. Disgusting. Mar 17, 2014 at 1:15

I was extremely tired when I wrote this.

I might have gone too far with projections idea, anyway, here's the projection the program uses. Basically like projecting earth onto a cube and then unfolding it. Besides, in this projection, the shadow is made of straight lines.
The program uses current date/time, and outputs a PPM file on stdout.

import Data.Time.Clock
import Data.Time.Calendar
import Control.Applicative
import Data.Fixed
import Data.Maybe

earth :: [[Int]]
earth = [[256],[256],[256],[256],[64,1,1,2,1,5,14,16,152],[56,19,3,27,1,6,50,1,2,1,90],[53,6,1,11,2,36,26,1,2,1,16,2,1,1,2,1,24,4,66],[47,2,5,14,4,35,22,7,54,2,1,3,60],[38,1,2,2,3,1,6,1,2,1,2,7,6,1,1,33,24,3,3,1,56,2,60],[34,2,1,4,2,1,3,1,1,3,3,2,15,3,3,29,57,5,19,1,2,11,17,1,1,1,34],[40,3,10,2,1,8,16,27,54,3,18,19,18,1,36],[33,6,5,3,2,3,1,3,2,2,1,5,16,21,1,2,53,2,10,1,6,19,1,7,4,3,9,2,33],[32,4,1,7,1,2,3,2,1,1,3,11,14,23,53,2,10,3,1,4,2,33,7,7,29],[8,5,25,10,5,3,2,14,10,2,1,18,1,2,31,6,18,1,7,4,1,60,22],[5,18,2,12,3,5,1,3,2,2,1,3,4,2,3,8,11,18,30,13,9,2,7,3,2,72,1,6,8],[4,36,2,1,1,4,3,7,1,4,3,9,8,15,34,18,2,2,2,17,1,78,4],[4,1,1,27,3,1,1,24,6,3,1,1,1,3,6,13,13,1,20,15,1,4,1,104,1],[3,31,1,24,1,2,4,8,10,9,12,6,18,7,3,7,1,1,2,99,3,2,2],[7,50,2,2,2,1,2,1,3,2,1,2,10,7,15,1,20,7,2,111,7,1],[4,35,1,15,9,1,1,3,4,1,12,5,34,8,3,110,10],[4,9,1,2,1,37,12,6,16,3,34,8,3,96,5,6,13],[6,6,1,1,8,32,12,6,3,1,49,9,4,2,1,86,1,3,4,2,19],[9,2,1,1,11,31,11,11,40,1,8,1,2,4,5,83,12,3,20],[8,1,16,33,9,11,39,2,8,1,2,3,3,83,13,5,19],[28,33,5,12,40,2,7,3,6,62,1,19,13,5,20],[27,36,2,15,34,3,2,2,6,71,1,22,11,2,22],[30,21,1,11,2,16,33,3,1,4,2,72,1,24,1,1,9,1,23],[31,21,1,26,39,4,1,98,1,1,33],[31,42,7,1,40,100,1,1,33],[33,25,2,15,4,4,35,102,36],[33,23,2,1,2,14,8,1,36,27,1,9,1,61,3,1,33],[33,26,5,14,42,10,1,11,2,2,2,7,3,5,1,9,1,44,38],[33,26,1,2,1,9,2,1,45,7,1,2,2,9,8,6,2,6,1,53,4,2,33],[33,26,1,4,1,6,44,8,6,2,3,7,9,5,3,56,1,1,4,3,33],[33,37,45,8,7,2,3,6,2,4,3,6,4,53,43],[33,36,46,6,6,1,4,1,2,2,3,16,3,47,1,5,8,2,34],[34,34,46,7,11,1,3,2,2,16,3,45,6,2,8,1,35],[34,33,48,5,11,1,4,1,4,16,2,49,3,2,6,2,35],[35,32,54,8,17,60,5,2,4,4,35],[36,30,50,12,18,60,8,2,1,1,38],[38,27,50,15,16,61,6,2,41],[38,25,51,18,3,4,6,62,6,1,42],[39,1,1,17,2,3,51,93,49],[40,1,1,11,9,2,49,31,1,10,2,50,49],[40,1,2,9,10,2,48,33,1,10,2,49,49],[41,1,2,8,11,1,47,34,2,10,5,44,50],[42,1,2,7,58,36,1,11,2,1,8,36,51],[46,6,58,36,2,15,7,34,2,1,49],[46,6,12,2,43,38,2,14,7,2,1,12,1,15,55],[46,6,5,2,7,2,41,38,2,14,10,10,4,10,59],[47,6,3,3,10,3,38,37,3,12,11,8,6,9,2,1,57],[49,10,51,38,3,9,13,7,8,9,9,2,48],[51,7,51,40,2,7,15,6,9,1,1,8,8,2,48],[55,7,47,41,1,6,17,4,12,8,8,1,49],[57,5,47,42,1,2,20,4,13,8,9,1,47],[59,3,8,1,38,43,22,4,13,1,2,4,10,2,46],[60,2,6,5,38,41,1,4,18,3,17,3,10,2,46],[61,2,1,1,2,3,1,7,34,45,18,2,18,1,60],[63,1,2,13,33,44,22,1,12,1,16,3,45],[66,14,33,43,22,1,13,1,14,1,1,1,46],[66,18,30,4,1,1,5,30,34,1,2,2,9,3,50],[66,19,43,27,34,2,2,1,7,3,52],[65,20,43,26,36,2,1,2,5,5,51],[65,21,42,24,39,3,4,7,2,1,1,1,1,1,44],[56,1,7,23,41,16,1,6,41,2,4,6,7,1,44],[64,25,39,16,1,5,42,3,4,5,2,1,8,1,2,1,37],[64,29,35,22,43,3,1,1,2,3,2,1,1,1,2,1,1,2,1,7,6,1,27],[63,31,35,20,45,2,11,1,9,7,4,2,26],[64,32,34,19,67,1,2,6,1,2,28],[65,31,34,12,1,6,48,4,18,6,31],[65,31,34,19,54,2,1,2,2,1,10,2,2,1,30],[66,29,36,14,1,3,57,1,19,2,28],[66,29,36,14,1,4,63,1,42],[67,27,36,15,1,4,63,5,3,2,33],[67,26,37,20,5,2,53,2,1,4,4,2,33],[68,25,37,20,4,3,52,9,3,3,32],[70,23,36,20,3,4,53,11,1,4,31],[71,22,37,17,5,4,51,18,31],[71,22,37,16,7,3,50,20,30],[71,21,39,15,6,3,5,1,42,24,29],[71,20,40,15,6,3,47,26,28],[71,17,43,15,6,3,46,28,27],[71,16,45,13,8,1,48,27,27],[71,16,45,12,58,28,26],[71,16,45,12,58,28,26],[70,16,47,10,59,28,26],[70,15,49,9,60,27,26],[70,14,50,7,62,7,6,13,27],[70,13,51,6,63,6,8,1,1,9,28],[70,10,138,10,28],[69,12,139,7,29],[69,11,141,5,19,3,8],[69,8,167,3,9],[69,8,166,1,1,1,10],[70,5,149,2,16,2,12],[69,6,166,3,12],[68,6,166,2,14],[68,5,166,3,14],[68,6,182],[67,6,183],[68,4,184],[68,4,6,2,176],[69,4,183],[70,5,20,1,160],[256],[256],[256],[256],[256],[256],[78,1,1,1,109,1,65],[75,2,115,1,23,1,39],[72,3,80,1,1,5,20,42,32],[74,1,70,1,4,21,5,52,2,1,25],[67,1,2,2,1,4,64,28,4,62,21],[69,9,34,1,1,1,1,1,1,1,2,48,3,69,15],[50,1,5,1,16,5,34,130,14],[32,1,1,2,4,1,3,1,4,29,32,128,18],[20,1,1,54,32,128,20],[17,49,34,137,19],[9,1,2,54,20,4,6,143,17],[16,51,18,5,10,135,21],[11,1,4,54,25,140,21],[12,66,4,155,19],[12,231,13],[0,6,9,5,2,234],[0,256],[0,256]]
main = do
mapM_ line [0..299]
where
putStrLn "P3"
putStrLn "# Some PPM readers expect a comment here"
putStrLn "400 300"
putStrLn "2"
line y = mapM_ (\x -> pixel x y >>= draw) [0..399]
where
draw (r, g, b) = putStrLn $(show r) ++ " " ++ (show g) ++ " " ++ (show b) pixel x y = fromMaybe (return (1, 1, 1))$
mapRegion (\x y -> (50, -x, y)) (x - 50) (y - 50)
<|> mapRegion (\x y -> (-x, -50, y)) (x - 150) (y - 50)
<|> mapRegion (\x y -> (-x, y, 50)) (x - 150) (y - 150)
<|> mapRegion (\x y -> (-50, y, -x)) (x - 250) (y - 150)
<|> mapRegion (\x y -> (y, 50, -x)) (x - 250) (y - 250)
<|> mapRegion (\x y -> (y, -x, -50)) (x - 350) (y - 250)
where
mapRegion f x y = if x >= -50 && y >= -50 && x < 50 && y < 50 then
Just $fmap (worldMap . shade) getCurrentTime else Nothing where t (x, y, z) = (atan2 y z) / pi p (x, y, z) = asin (x / (sqrt$ x*x+y*y+z*z)) / pi * 2
rotate o (x, y, z) = (x, y * cos o + z * sin o, z * cos o - y * sin o)
tilt o (x, y, z) = (x * cos o - y * sin o, x * sin o + y * cos o, z)
shade c = ((t $rotate yearAngle$ tilt 0.366 $rotate (dayAngle - yearAngle)$ f x y)) mod' 2 > 1
where
hour=(realToFrac$utctDayTime now)/(60*60) mnlong=280.460+0.9856474*time mnanom=(357.528+0.9856003*time)*d eclong=(mnlong+1.915*sin(mnanom)+0.020*sin(2*mnanom))*d oblqec=(23.439-0.0000004*time)*d ra=let num=cos(oblqec)*sin(eclong) den=cos(eclong) in if den<0 then atan(num/den)+pi else atan(num/den) dec=asin(sin(oblqec)*sin(eclong)) gmst =6.697375+0.0657098242*time+hour lmst=(gmst*15*d)+long ha=(lmst-ra) el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha)) in el<=0 td x = fromIntegral x :: Double keep="NSEW"++['0'..'9'] pixel p dk=if dk && pnotElemkeep then if p==' ' then '#' else '%' else p showMap t= do let w=length(worldmap!!0) h=length worldmap putStrLn (worldmap!!0) putStrLn (worldmap!!1) mapM_(\y->do mapM_(\x->let lat=(0.5-td y/td h)*pi long=(0.5-td x/td w)*tau in putStr [pixel ((worldmap!!(y+2))!!x) (dark lat long t)]) [0..(w-1)] putStrLn "") [0..(h-4)] putStrLn (last worldmap) main = do {t<-getCurrentTime; showMap t} worldmap=[ "180 150W 120W 90W 60W 30W 000 30E 60E 90E 120E 150E 180", "| | | | | | | | | | | | |", "+90N-+-----+-----+-----+-----+----+-----+-----+-----+-----+-----+-----+", "| . _..::__: ,-\"-\"._ |7 , _,.__ |", "| _.___ _ _<_>!(._.-. / _._ _ ,_/ ' '-._.---.-.__|", "|.{ \" \" -==,',._\\{ \\ / {) / _ \">_,-'  mt-2_|", "+ \\_.:--. ._ )^-. \"' , [_/( __,/-' +", "|'\"' \\ \" _L oD_,--' ) /. (| |", "| | ,' _)_.\\\\._<> 6 _,' / ' |", "| . / [_/_' \"( <'} ) |", "+30N \\\\ .-. ) / -'\"..' :._ _) ' +", "|  \\ ( ( / :\\ > \\ ,-^. /' ' |", "| ._, \"\" | \\' \\| ?_) {\\ |", "| =.---. ._._ ,' \" |' ,- '. |", "+000 | -._ | / :<_|h--._ +", "| ( > . | , =.__.-'\\ |", "| . / | |{| ,-.,\\ .|", "| | ,' \\ / ' ,\" \\ |", "+30S | / |_' | __ / +", "| | | '-' -' \\.|", "| |/ \" / |", "| \\. ' |", "+60S +", "| ,/ ______._.--._ _..---.---------._ |", "| ,-----\"-..?----_/ ) _,-'\" \" ( |", "|.._( -----' -|", "+90S-+-----+-----+-----+-----+----+-----+-----+-----+-----+-----+-----+", "Map 1998 Matthew Thomas. Freely usable as long as this line is included"]  Example output, from a more interesting time of year (we're near the equinox, so Wander Nauta's rectangular blobs are fairly accurate :) ) - this is for Jan 16 13:55:51 UTC 2014: 180 150W 120W 90W 60W 30W 000 30E 60E 90E 120E 150E 180 | | | | | | | | | | | | | %90N%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %##########%#%%%%%%%%##%%%%%%%#######%7#######%#####%%%%%#############% %##%%%%%#%#%%%%%%%%%%%%%%####%########%%%#####%%#%%%##%##%%%%%%%%%%%%%% %%%#####%#%#%%%%%%%%%%%##%##%#%%#####%#%#%%%%%%#%################%%%2%% %#%%%%%%%#######%%%#%%%%%#%%######, [_/( ##############%%%%%%#% %%%%#####%#########%####%%##### oD_,--' ####%#####%%#%%###% %#########%###########%%##### _)_.\\._<> 6 ######%%%#%##%###% %#########%%#########%###### [_/_' "( ###%%%##%######% %30N#######%%####%%%#%##### / -'"..' :._ ###%%##%#######% %###%########%##%##%%##### / :\ > \ ,-^. #%%#%#########% %#############%%%%###%%### | \' \| ?_)##%%#########% %################%%%%%%%# ._._ ,' " |' %%#%%########% %000###############%####-._ | / :<_%%%%%%######% %##################%#### > . | , =.%%%%%%%#####% %###################%%# / | |{| %%%%%#####%% %####################%# ,' \ / ' ,"#####%#####% %30S#################% / |_' | %%##%#####% %####################% | '-'##%%%###%%% %####################|/ ##%####%#% %####################\. #####%##% %60S################ ########% %################## ,/ ______._.--._ _..---.-------%%%%###% %####%%%%%%%%%%%%%--_/ ) _,-'" " ##%##% %%%%%########### -----' ##%%% %90S%%%%%%%%%----+-----+-----+----+-----+-----+-----+-----+-----+----%% Map 1998 Matthew Thomas. Freely usable as long as this line is included  [1] they're the same as you'll find elsewhere, except without the extra work to keep degrees between 0 and 360, hours between 0 and 24, and radians between 0 and 2pi. I think those are holdovers from the days we used slide rules; trig functions work just fine outside those ranges... • Brilliant! I love that you can still see the map through the 'dark'. Also, the math looks solid. Could you add the date you used for the example, so other people can compare their solutions to yours? Mar 17, 2014 at 0:24 • I see that you've added the date, thanks! Mar 17, 2014 at 0:26 • Yep. I like that in the image you can clearly see that it's winter in the Northern Hemisphere, made it easier to believe that date! I'm glad you posted an answer before me, saved me endless faff trying to do a golfed version, there's no way I'll beat yours for that. Mar 17, 2014 at 0:29 . Bash, 882* characters This is my second entry, this time in the Aesthetics, Weird tech, Fun and Short code categories. It's inspired by Ram Narasimhan's entry and Peter Taylor's comment. The script first generates a low-res texture of the world, bundled as base64-encoded data. It then generates 24 PovRay-scenes containing a sphere with that texture, each one rotated to 'face the sun'. Finally, the frames are combined into a GIF animation using ImageMagick. This means you'll have to have both PovRay and ImageMagick installed for the script to work - feel free to ignore this entry if you think that should disqualify it. Like Ram's entry, and my first entry, this does not account for seasonal change, which means it's not very precise. It is, however, shorter, prettier and more precise than my first entry - and unlike Ram's entry, the map data and the code for generating the GIF animation are included.  echo ' iVBO Rw0KGgoAAAA NS UhE U g AAAEgAAAA kAQMAAAAQFe4lAAAABlB MVEUAFFwAbxKgAD63 AAAA AWJLR0 QAiAUdSAAAAAlwSFlzAAALEwAACx MB AJqcGAAAAAd0SU1FB9 4DE hUWI op Fp5MAAADDSURBVBhXrcYhTsNQGADgr3ShE4Qi h4BeYQFBgqAJN8Lh +r jBb rArIJHPobgAgkzgeSQkVHT7MWThAHzq44 /j/jezy6jSH M6fB gd 9T Nbxdl99R4Q+XpdNRISj4dlFRCz oI11FxIpup4uIRDe5 fokp0Y2W25jQFDfrGNGsDNsoqBaGj34D2 bA7TcAwnmRoDZM 5tLkePUJb6uIT2rEq7hKaUhUHCXWpv7Q PqEv1rsuoc7X RbV Bn2d kGTYKMQ3C7H8z2+wc/eMd S QW39v8kAAA AA SUVOR K5CYII='|base64 \ -di>t;for X in {0..23};do R=$((90-(\
$X*15) )); echo "camera{location <0, 0, -5> angle 38 } light_source{ <0,0, -1000> rgb < 2,2, 2>} sphere {<0 ,0,0> 1 pigment { /**/ image_map{\"t\" map_type 1}} rotate <0,$R,0>
}">s               ;povray             +Is +H300\
+Of$X.png +W400 mogrify -fill white \ -annotate +0+10 "$X:00" \
-gravity           south                    f$X.png done; convert -delay \ 100 -loop 0$(ls f*  \
|sort               -V)                  ani.gif
exit;


As a bonus, here's a GIF that uses NASA's Blue Marble image instead of the space-saving 1-bit texture, i.e. what the result would have looked like without any size restriction: http://i.imgur.com/AnahEIu.gif

*: 882 characters not counting decorative whitespace, 1872 characters total.

• +1 for making it all self-contained. And also for creating the self-referential code which itself looks like a world map. Nice work. Mar 19, 2014 at 0:49
• argh! multiple entries. Now I'm losing my excuse for not doing a mad one... Mar 19, 2014 at 1:29
• That's sick! Love it. Mar 19, 2014 at 8:04
• Heh, I realize now that I could've used the formatting of the code itself as the source for the world map (with spaces = ocean, everything else = land) and actually get better resolution with less characters. Oh well... Mar 19, 2014 at 11:21
• This looks like Minecraft. Oct 20, 2014 at 20:57

I decided to kick off the contest with an entry of my own, in the short code category. It's 923 characters long, not counting newlines.

C: 923 characters

Here's the code:

i;j;w=160;f=40;t;b;p;s;e;k;d=86400;q=599;
char* m="M('+z EDz :!#\"!*!8S$[\"!$!#\"\")\"!3R)V$'!!()1M./!F)\"!!!!)'/GE5@\"\"!&%.3&,Y$D\"!!%$)5i\"\"\"F\"%&&6%!e'A#!#!!#&$5&!f&A'$*\"5&!c-#'3''8\"$!!#\"U'\"=5$'8#$$\"S(#=7!*5\"!\"#['!A@6#!^H=!#6bH;!!!\"6_!!I;<&!&\"!!\"F\"!I8;&\"#\"&#\"C#\"I7<%#!\"/\"BP5=*,\"=#\"!L4A%&\"\"G\"\"\"#M1@)*F\"%P/@,!N#!S(E;!@W'E=!!!<Y&D7!&!\"7\\D8!)4_C8!('&#&!!a&@9!&(%&g>9!*#(%h\">:!!-\"(%&!b!&5:!\"+\"(!!#!!!c+5<-!'!'!#!e)5:.!(!&!\"\"e,:25!!!\"!\"\"h-;07#\"h.9/:\"\"!!#\"a17-;'!\"!!\"!X46,<\"%\"&\\45,>#&!$$#!W45,C!!!'!\"!$!V26,H\"#!$!\"!\"!S17-#!A!!#\"!_07,\"#A&!\".7+#\"A*.!Q.7*$\">/^-9)$\"=0^*<)$!>1]*<(D1])>&E2\\)>&F&!)\\)@#G$%(\\'w%]'x#,\"P%z .\"P%z .!R$z -\"S$z b#z c#z d#z 3"; main(){ t=(time(0)%d*160)/d; printf("P2\n%d 62\n5\n",w); for(;i<q;i++){ for(j=m[i]-' ';j>0;j--){ p=k%w,s=(t-f),e=(t+f); printf("%c ","1324"[b*2+((p>s&&p<e)||(p>s+w&&p<e+w)||(p>s-w&&p<e-w))]); k++; } b=!b; } }  Here's how it works: A crude bitmap of the world* is run-length encoded as a string. Every character in the string represents a run of either land or sea pixels. Long runs of sea are split into a run of sea, then 0 land pixels, then another run of sea, to avoid including unprintable characters in the string. The Python script I wrote to convert PBM files into this format is here. I then use time() to find out how many seconds have passed in Greenwich since midnight, 1 January 1970. I modulo that to find out how many seconds have passed there today, using that information to position the light portion of the map more-or-less accordingly (I hope). Correctness is a joke. There's no math at all. The code assumes the earth is a cylinder (block-shaped day/night), that the sun is directly above the equator (no summer/winter), and that you like the color gray (no color). On the plus side, I do draw continents. The output is in Portable Graymap (PGM) format, which can then be converted to PNG by something like ImageMagick or the GIMP. Here's an example output, converted to PNG (larger version): *: The entire world except Antarctica, but who lives there anyway... • Nice output, it would be better if the dark area was curved though – qwr Mar 17, 2014 at 0:46 • Yep! Having the output be curved means you would have to do some trigonometry though, which would make it much longer. (Or I guess you could just round the correct corners to make it look right-ish...) Mar 17, 2014 at 0:50 • @WanderNauta There'll be some severely pissed off penguins that will feel cheated that you didn't take them into consideration, bro... Mar 17, 2014 at 10:43 • @WallyWest If you live on the poles, like your penguins, you don't need this program - midnight sun and all. Mar 17, 2014 at 10:48 Haskell - it's Hammer time. I did another one. Adapted from my previous version, this one uses an oblique Hammer projection to show both poles at the same time (in fact you're seeing the whole earth in each frame). Just for added weirdness, instead of using a bitmap directly, I sampled the earth along a spiral to give approximately equal area coverage; this is what allows me to distort the earth and rotate it easily. Hammer's projection is equal area too; my idea was that pairing these two things would lead to less distortion when I fill in the gaps. I display a graticule on the projection too, using Bresenham's algorithm to draw the lines. The globe and terminator line both move over the course of the day. Edited to change the image to a higher resolution but (deliberately) coarser underlying bitmap, so you can see the effect of the spiral. This uses 5000 points (sampled from ~260000), equivalent to a 50x100 bitmap, but giving more resolution to the equator than to the poles. To use the code, compile with ghc, run with an optional numeric parameter which is the hour offset; files are generated like 'earth0.pgm', 'earth1.pgm'. import System.Environment import Data.List (intercalate,unfoldr) import qualified Data.Set as Set import Data.List.Split import Data.List import Data.Maybe (catMaybes) import qualified Data.Map as Map import Data.Time import Debug.Trace d=pi/180 tau=2*pi m0=UTCTime(fromGregorian 2000 1 1)(secondsToDiffTime(12*60*60)) dark::Double->Double->UTCTime->Bool dark lat long now = let time=(realToFrac$diffUTCTime now m0)/(60*60*24)
hour=(realToFrac$utctDayTime now)/(60*60) mnlong=280.460+0.9856474*time mnanom=(357.528+0.9856003*time)*d eclong=(mnlong+1.915*sin(mnanom)+0.020*sin(2*mnanom))*d oblqec=(23.439-0.0000004*time)*d ra=let num=cos(oblqec)*sin(eclong) den=cos(eclong) in if den<0 then atan(num/den)+pi else atan(num/den) dec=asin(sin(oblqec)*sin(eclong)) gmst =6.697375+0.0657098242*time+hour lmst=(gmst*15*d)+long ha=(lmst-ra) el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha)) in el<=0 infill(open, known)= if null open then known else infill gen where neighbours (x,y)=catMaybes$ map ((flip Map.lookup) known) [(x+1,y),(x-1,y),(x,y+1),(x,y-1),(x+1,y+1),(x-1,y+1),(x-1,y-1),(x-1,y-1)]
vote a= if null a then Nothing
else Just ((sum a)div(length a))
fill x (open',  known')=
case vote (neighbours x) of
Nothing->(x:open',known')
Just c->(open',(x,c):known')
gen=(\(o,k)->(o,Map.fromList k))$foldr fill ([], Map.toList known) open mpoint (a,b)=case a of Nothing->Nothing;Just c->Just(c,b) grid w h n g lut= map (\y->map (\x->if Set.member (x,y) g then 3 else case Map.lookup (x,y) lut of Nothing->7;Just c->c) [1..w]) [1..h] unknowns w h lut=concatMap (\y->concatMap (\x->let z=1-(2*x//w-1)^2-(2*y//h-1)^2 in case Map.lookup (x,y) lut of Nothing->if z<0 then [] else [(x,y)];_->[]) [1..w]) [1..h] main=do args <- getArgs let off = if null args then 0 else read(args!!0) actual <- getCurrentTime let now=((fromIntegral off)*60*60) addUTCTime actual let tod=realToFrac(utctDayTime now)/86400+0.4 let s=5000 let w=800 let h=400 let n=6 -- pbm <- readFile "earth.pbm" -- let bits=ungrid s$parsepbm pbm
let bits=[0,23,4,9,1,3,1,2,6,10,1,10,4,1,3,7,10,7,4,2,2,1,2,6,12,1,1,2,1,5,4,1,8,1,3,
1,21,7,2,2,35,1,4,3,2,2,2,2,16,1,25,1,2,8,1,4,1,2,13,3,2,1,26,1,1,10,3,3,8,
2,3,6,1,3,25,2,1,10,15,5,1,6,2,3,30,10,15,19,32,11,16,20,35,11,1,2,14,22,27,
1,8,14,16,22,2,1,22,1,1,2,1,1,2,1,2,1,3,16,14,25,1,2,21,1,6,1,2,1,1,2,3,17,
14,26,1,2,1,1,26,1,1,3,3,1,1,19,13,28,4,1,26,6,6,21,11,35,40,21,11,37,41,20,
2,4,4,1,1,39,19,1,6,1,16,19,2,4,5,40,18,2,7,1,17,19,1,1,1,1,1,2,3,46,7,1,5,
4,25,16,3,1,1,3,5,44,1,4,5,4,3,6,4,1,19,22,5,46,2,3,4,6,2,9,22,22,2,50,1,5,
2,1,1,6,1,8,24,15,5,1,2,51,2,5,1,1,1,5,1,10,23,14,9,55,1,4,2,17,16,1,4,14,9,
57,4,1,3,17,13,20,11,54,2,1,3,1,2,20,12,18,13,47,4,3,8,21,10,17,15,44,5,1,1,
4,1,3,2,22,10,15,16,46,4,3,1,2,2,25,9,17,15,47,1,1,3,30,9,18,13,46,2,1,4,25,
2,1,11,16,13,46,8,24,2,2,9,16,11,45,12,22,1,3,7,17,10,45,12,21,1,3,7,19,8,
43,12,25,6,19,8,41,12,25,5,20,7,40,11,25,4,20,6,40,5,3,2,48,6,38,3,54,4,30,
1,6,2,55,2,29,1,5,1,53,3,28,1,55,3,49,1,30,2,76,1,284,3,4,1,15,1,17,10,1,9,
7,1,13,21,4,4,1,2,6,17,2,8,3,63]
let t(phi,lambda)=unitsphere$rx (-pi/4)$rz (-tod*4*pi)$sphereunit(phi, lambda) let hmr=(fmap (\(x,y)->(floor((fl w)*(x+4)/8),floor((fl h)*(y+2)/4)))).hammer.t let g=graticule hmr n let lut = Map.fromList$ catMaybes $map mpoint$map (\((lat,long),bit)->(hmr(lat,long),bit*4+2-if dark lat long now then 2 else 0))  $zip (spiral s) (rld bits) -- let lut = Map.fromList$ catMaybes $map mpoint$map (\((lat,long),bit)->(hmr(lat,long),bit))$zip (spiral s) (rld bits) let lut' = infill ((unknowns w h lut), lut) let pgm = "P2\n"++((show w)++" "++(show h)++" 7\n")++(intercalate "\n"$ map (intercalate " ")$chunksOf 35$ map show(concat$grid w h n g lut'))++"\n" writeFile ("earth"++(show off)++".pgm") pgm fl=fromIntegral spiral::Int->[(Double,Double)] spiral n=map (\k-> let phi=acos(((2*(fl k))-1)/(fl n)-1) in rerange(pi/2-phi,sqrt((fl n)*pi)*phi)) [1..n] rld::[Int]->[Int] rld bits=concat$rld' (head bits) (tail bits)
where
rld' bit []=[]
rld' bit (run:xs) = (replicate run bit):(rld' (case bit of 1->0;_->1) xs)
rle::[Int]->[Int]
rle bits=(head bits):(map length$group bits) sample::Int->Int->Int->[(Int,Int)] sample n w h = map (\(phi, theta)->((floor((fl h)*((phi-(pi/2))/pi)))modh, (floor((fl w)*(theta-pi)/(tau)))modw ))$ spiral n
ungrid::Int->[[Int]]->[Int]
ungrid n g = rle $map (\(y, x)->(g!!y)!!x) (sample n w h) where w = length$head g
h = length g
parsepbm::[Char]->[[Int]]
parsepbm pbm=
[width, height] = map read$words (head$drop 1 header)
d = ((map read).concat.(map words)) rest
in chunksOf width d
rerange(phi,lambda)
| abs(phi)>pi = rerange(phi - signum(phi)*tau, lambda)
| abs(phi)>pi/2 = rerange(phi-signum(phi)*pi, lambda+pi)
| abs(lambda)>pi = rerange(phi, lambda - signum(lambda)*tau)
| otherwise = (phi, lambda)
laea(phi,lambda)=if isInfinite(z) then Nothing else Just (z*cos(phi)*sin(lambda),z*sin(phi)) where z=4/sqrt(1+cos(phi)*cos(lambda))
hammer(phi,lambda)=case laea(phi, lambda/2) of Nothing->Nothing; Just(x,y)->Just (x, y/2)
bresenham :: (Int, Int)->(Int, Int)->[(Int, Int)]
bresenham p0@(x0,y0) p1@(x1,y1)
| abs(dx)>50||abs(dy)>50=[]
| x0>x1 = map h$bresenham (h p0) (h p1) | y0>y1 = map v$ bresenham (v p0) (v p1)
| (x1-x0) < (y1-y0) = map f$bresenham (f p0) (f p1) | otherwise = unfoldr (\(x,y,d)->if x>x1 then Nothing else Just((x,y),(if 2*(d+dy)<dx then(x+1,y,d+dy)else(x+1,y+1,d+dy-dx)))) (x0,y0,0) where h(x,y)=(-x,y) v(x,y)=(x,-y) f(x,y)=(y,x) dx=x1-x0 dy=y1-y0 globe n k= (concatMap (\m->map (meridian m) [k*(1-n)..k*(n-1)]) [k*(1-2*n),k*(2-2*n)..k*2*n]) ++(concatMap (\p->map (parallel p) [k*(-2*n)..k*2*n]) [k*(1-n),k*(2-n)..k*(n-1)]) where meridian m p=(radians(p,m),radians(p+1,m)) parallel p m=(radians(p,m),radians(p,m+1)) radians(p,m)=rerange((p//(k*n))*pi/2,(m//(k*n))*pi/2) graticule f n=Set.fromList$ concatMap (\(a,b)->case (f a,f b) of (Nothing,_)->[];(_,Nothing)->[];(Just c,Just d)->bresenham c d) (globe n 4)
rx theta (x,y,z) = (x, y*(cos theta)-z*(sin theta), y*(sin theta)+z*(cos theta))
ry theta (x,y,z) = (z*(sin theta)+x*(cos theta), y, z*(cos theta)-x*(sin theta))
rz theta (x,y,z) = (x*(cos theta)-y*(sin theta), x*(sin theta)+y*(cos theta), z)
sphereunit (phi, theta) = (rz theta (ry (-phi) (1,0,0)))
unitsphere (x,y,z) = (asin z, atan2 y x)
x//y=(fromIntegral x)/(fromIntegral y)

• This is... madness. Absolute madness. I love it. Mar 27, 2014 at 20:14
• It tooks me so long to get the transform right on this one. The spiral bitmap alone is hilarious. Mar 27, 2014 at 20:16
• re "it does scale up quite well" - do you mean you could easily make a higher-resolution GIF? Mar 27, 2014 at 20:23
• Yes. quite a lot larger. I took that code out for a bit while I messed with the graticule, when I get time a bit later I'll put it back in and upload a better image Mar 27, 2014 at 20:25
• there you go - nice big image of crazytown Mar 28, 2014 at 2:27

C, using pnm images

Late answer, focusing on correctness and aesthetics. The output is a blend of two input images (day.pnm and night.pnm), including a stripe of twilight. I'm using images based on NASAs blue marble here.

The code uses my own img.h for clarity (just imagine it being included verbatim in the .c for strict rule compliance...). Everything in there is implemented via C macros. The animations are built with imagemagicks convert from multiple frames - the program itself will only output static images. Code is below.

Now: (Aug 13, ~13:00 CEST)

One day: (Jan 1st)

One year: (12:00 UTC)

sun.c

  #include <math.h>
#include <time.h>

#include "img.h"

#ifndef M_PI
#define M_PI 3.14159265359
#endif

double deg2rad(double x) {return x / 180.0 * M_PI;}
double rad2deg(double x) {return x * 180.0 / M_PI;}

double elevation(double latitude, double longitude, int yday, int hour, int min, int sec)
{
double fd = (hour + (min + sec / 60.0) / 60.0) / 24.0;
double fyd = 360.0 * (yday + fd) / 366.0;

double m = fyd - 3.943;
double ta = -1.914 * sind(m) + 2.468 * sind(2 * m + 205.6);
double hourangle = (fd - 0.5) * 360.0 + longitude + ta;
double decl = 0.396 - 22.913 * cosd(fyd) + 4.025 * sind(fyd) - 0.387 * cosd(2 * fyd) + 0.052 * sind(2 * fyd) - 0.155 * cosd(3 * fyd) + 0.085 * sind(3 * fyd);

return asind(cosd(hourangle) * cosd(decl) * cosd(latitude) + sind(decl) * sind(latitude));
}

int main(int argc, char* argv[])
{
Image day, night, out;
int x, y;
time_t t = time(0);
struct tm* utc = gmtime(&t);
int yday = utc->tm_yday, hour = utc->tm_hour, min = utc->tm_min, sec = utc->tm_sec;

for(y = 0; y < day.height; ++y)
{
double latitude = 90.0 - 180.0 * (y + 0.5) / day.height;
for(x = 0; x < day.width; ++x)
{
double longitude = -180.0 + 360.0 * (x + 0.5) / day.width;
double elev = elevation(latitude, longitude, yday, hour, min, sec);
double nf = elev > -0.8 ? 0.0 : elev > -6.0 ? 0.5 : 1.0;
double df = 1.0 - nf;
Color dc = imgGetColor(day, x, y);
Color nc = imgGetColor(night, x, y);
imgDotC3(out, x, y, df * dc.r + nf * nc.r, df * dc.g + nf * nc.g, df * dc.b + nf * nc.b);
}
}
imgSave(out, "out.pnm");
}


img.h

  #include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef struct
{
unsigned char r;
unsigned char g;
unsigned char b;
} Color;

typedef struct
{
Color* data;
int width;
int height;
Color c;
} Image;

#define imgCreate(img, w, h)           {\
int length;\
(img).width = (w);\
(img).height = (h);\
length = (img).width * (img).height * sizeof(Color);\
(img).data = malloc(length);\
memset((img).data, 0, length);\
(img).c.r = (img).c.g = (img).c.b = 0;\
}

#define imgDestroy(img)                {\
free((img).data);\
(img).width = 0;\
(img).height = 0;\
(img).c.r = (img).c.g = (img).c.b = 0;\
}

#define imgSetColor(img, ur, ug, ub)   {\
(img).c.r = (ur);\
(img).c.g = (ug);\
(img).c.b = (ub);\
}

#define imgDot(img, x, y)              {\
(img).data[(int)(x) + (int)(y) * (img).width] = (img).c;\
}

#define imgDotC3(img, x, y, ur, ug, ub) {\
(img).data[(int)(x) + (int)(y) * (img).width].r = (ur);\
(img).data[(int)(x) + (int)(y) * (img).width].g = (ug);\
(img).data[(int)(x) + (int)(y) * (img).width].b = (ub);\
}

#define imgDotC(img, x, y, c)          {\
(img).data[(int)(x) + (int)(y) * (img).width] = (c);\
}

#define imgGetColor(img, x, y)         ((img).data[(int)(x) + (int)(y) * (img).width])

#define imgLine(img, x, y, xx, yy)     {\
int x0 = (x), y0 = (y), x1 = (xx), y1 = (yy);\
int dx =  abs(x1 - x0), sx = x0 < x1 ? 1 : -1;\
int dy = -abs(y1 - y0), sy = y0 < y1 ? 1 : -1;\
int err = dx + dy, e2;\
\
for(;;)\
{\
imgDot((img), x0, y0);\
if (x0 == x1 && y0 == y1) break;\
e2 = 2 * err;\
if (e2 >= dy) {err += dy; x0 += sx;}\
if (e2 <= dx) {err += dx; y0 += sy;}\
}\
}

#define imgSave(img, fname)            {\
FILE* f = fopen((fname), "wb");\
fprintf(f, "P6 %d %d 255\n", (img).width, (img).height);\
fwrite((img).data, sizeof(Color), (img).width * (img).height, f);\
fclose(f);\
}

FILE* f = fopen((fname), "rb");\
char buffer[16];\
int index = 0;\
int field = 0;\
int isP5 = 0;\
unsigned char c = ' ';\
while(field < 4)\
{\
do\
{\
if(c == '#') while(c = fgetc(f), c != '\n');\
} while(c = fgetc(f), isspace(c) || c == '#');\
index = 0;\
do\
{\
buffer[index++] = c;\
} while(c = fgetc(f), !isspace(c) && c != '#' && index < 16);\
buffer[index] = 0;\
switch(field)\
{\
case 0:\
if (strcmp(buffer, "P5") == 0) isP5 = 1;\
else if (strcmp(buffer, "P6") == 0) isP5 = 0;\
else fprintf(stderr, "image format \"%s\" unsupported (not P5 or P6)\n", buffer), exit(1);\
break;\
case 1:\
(img).width = atoi(buffer);\
break;\
case 2:\
(img).height = atoi(buffer);\
break;\
case 3:\
index = atoi(buffer);\
if (index != 255) fprintf(stderr, "image format unsupported (not 255 values per channel)\n"), exit(1);\
break;\
}\
field++;\
}\
imgCreate((img), (img).width, (img).height);\
if (isP5)\
{\
int length = (img).width * (img).height;\
for(index = 0; index < length; ++index)\
{\
(img).data[index].r = (img).data[index].g = (img).data[index].b = fgetc(f);\
}\
}\
else\
{\
fread((img).data, sizeof(Color), (img).width * (img).height, f);\
}\
fclose(f);\
}

• Very pretty! I like it. Aug 13, 2014 at 13:00

R: Using ggplot2 and Map projection

Inspired by @mniip's post, I decided to try using R's mapproj package, wherein we can orient the globe by specifying where the North Pole should be when computing the projection.

Based on current GMT time, I compute the longitude where is it currently noon and make that point the map's center. We are looking at Earth from the "Sun's point of view" so whatever is visible is in daylight.

Much of the code is just aesthetics. The only part I had to figure out was to compute the "noon Longitude," that is the longitude value where it was noon given any GMT time.

library(ggplot2);library(maps);library(ggmap)
world <- map_data("world")# a lat-long dataframe from the maps package
worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) +
geom_path(color="orange") +
theme(panel.background= element_rect("black"),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.grid.major = element_line(colour="blue", size=0.75),
panel.grid.minor = element_line(colour="blue")
)

#Create a function that takes in the current GMT time
print_3d_coordmap <- function (current_gmt_time) {
curr_gmt_mins <- as.POSIXlt(current_gmt_time)$hour*60 + as.POSIXlt(current_gmt_time)$min
noon_longitude <- 180 - (curr_gmt_mins * 360/1440)
#centered at wherever longitude where it is Noon now on (lat:equator)
worldmap + coord_map("ortho", orientation=c(0, noon_longitude, 0))
}

#test it out
print_3d_coordmap(Sys.time() + 7*60*60) # my location is 7 hours behind UTC


I then used the R animation package to generate 24 images and stitched them into one GIF.

• Looks great! Does it correctly handle summer and winter, though? I don't know R that well, but it looks like your equator is always at the center of the image. Mar 18, 2014 at 9:33
• Yes, you are correct. I did a quick and dirty implementation using the equator as the only latitude. (Focused on the animation instead.) The project has lots of features that I am not using. If there is a reference that shows how to vary latitudes with seasons, I will be happy to try it out. Mar 18, 2014 at 19:06
• Here's a reference in R, Ram - in fact this is where I translated the formulas in my entry from stackoverflow.com/questions/8708048/… Mar 18, 2014 at 20:15
• @bazzargh Thanks! Looks like I have a lot of learning about azimuth angles to do. Mar 19, 2014 at 0:46

JavaScript – by Martin Kleppe (http://aem1k.com/)

I want to stress that this is not my work, but the work of Martin Kleppe. I just think it fits so perfectly that it should not be missing here:

Online Demo (or just paste it into the console)

eval(z='p="<"+"pre>"/*        ######## */;for(y in n="zw24l6k\
4e3t4jnt4qj24xh2 x/*    *############### */42kty24wrt413n243n\
9h243pdxt41csb yz/*  #################### */43iyb6k43pk7243nm\
r24".split(4)){/*     *#################*   */for(a in t=pars\
eInt(n[y],36)+/*          ###############*   */(e=x=r=[]))for\
(r=!r,i=0;t[a/*               ############*   */]>i;i+=.05)wi\
th(Math)x-= /*                #############    */.05,0<cos(o=\
new Date/1e3/*                #########*       */-x/PI)&&(e[~\
~(32*sin(o)*/*                     ####*       */sin(.5+y/7))\
+60] =-~ r);/*                         *###    */for(x=0;122>\
x;)p+="   *#"/*                        #####  */[e[x++]+e[x++\
]]||(S=("eval"/*                      *##### */+"(z=\'"+z.spl\
it(B = "\\\\")./*      ###*           ####  */join(B+B).split\
(Q="\'").join(B+Q/*                  ###* */)+Q+")//m1k")[x/2\
+61*y-1]).fontcolor/*               ##   */(/\\w/.test(S)&&"#\
03B");document.body.innerHTML=p+=B+"\\n"}setTimeout(z)')//m1k\

• If it's not your work, you should turn your answer into a community wiki. Aug 13, 2014 at 14:28
• Thanks for the tip, never noticed that checkbox before. Done! Aug 13, 2014 at 14:31