Simulating Sunlight

Version 1: December 1999
Version 2: April 2010

Matthew P. Szudzik
Wolfram Research, Inc.
(mszudzik@wolfram.com)

Definitions

ToSequence[x___] := x

basis[m_] :=
Module[{a},
a =First[m];
DeleteCases[RowReduce[Map[a-#&,Rest[m]]],x_ /; x.x==0]
]

Options[Sunlight] = {Beam->False,ClipDistance->.0001,Orientation->0 Degree,Epilog->{},Prolog->{},RoomDimensions->{1,1},WindowPoints->{{0,0,0},{0,0,1},{0,1,1},{0,1,0}}};

Sunlight[date_,opts___] :=
Module[{beam,clip,dir,epilog,floor,len,orient,points,prolog,vector,wid,x,y,z},
beam = Beam /. {opts} /. Options[Sunlight];
clip = ClipDistance /. {opts} /.Options[Sunlight];
epilog = Epilog /. {opts} /.Options[Sunlight];
orient =Orientation /. {opts} /.Options[Sunlight];
dir =Direction /. {opts} /. Options[Sunlight];
If[dir=!=Direction,orient=dir+180 Degree];
points =Map[If[Length[#]===2,Prepend[#,0.],#]&,N[WindowPoints /. {opts} /.Options[Sunlight]],{-2}];
prolog = Prolog /. {opts} /.Options[Sunlight];
{x,y} = RoomDimensions /. {opts} /.Options[Sunlight];
len =RoomLength /. {opts};
wid =RoomWidth /. {opts};
If[len=!=RoomLength,x=len];
If[wid=!=RoomWidth,y=wid];
z =Max[Flatten[Map[#[[3]]&, points,{-2}]]];
vector=Function[{az,alt},N[-{Cos[orient-az]*Cos[alt],Sin[orient -az]*Cos[alt],Sin[alt]}]]@@(AstronomicalData["Sun",{#,date}] Degree&/@{"Azimuth","Altitude"});
Graphics3D[{ToSequence @@prolog,{SurfaceColor[RGBColor[.75,.85,1]],Map[Polygon,points,{-3}]},If[ vector[[3]]≥0 || (points=Flatten[{Map[If[(((#1 . #2)*#1) .vector)&[Cross @@Take[basis[#],2],{x/2,y/2,0}-#[[1]]] > 0,#,ToSequence[]]&,points,{-3}]},1])==={},ToSequence[],
floor =Map[#+vector*(0-#[[3]])/vector[[3]]&,points,{-2}];ToSequence[If[beam =!= True,ToSequence[], {SurfaceColor[RGBColor[1,1,.5]],MapThread[Polygon[{#1[[1]],#1[[2]],#2[[2]],#2[[1]]}]&,If[Depth[points]===4, {Flatten[#[[1]],1],Flatten[#[[2]],1]},#]&[Map[Partition[#,2,1,{1,1}]&,{points,floor},{-3}]]]}], {SurfaceColor[RGBColor[1,1,0]],Map[ Polygon,floor,{-3}],Map[If[vector[[1]] ≤ 0,ToSequence[],Polygon[#]]&,Map[#+vector*(x-#[[1]])/vector[[1]]&,points,{-2}],{-3}],
Map[If[vector[[1]] ≥ 0,ToSequence[],Polygon[#]]&,Map[#+vector*(0-#[[1]])/vector[[1]]&,points,{-2}],{-3}],
Map[If[vector[[2]] ≤ 0,ToSequence[],Polygon[#]]&,Map[#+vector*(y-#[[2]])/vector[[2]]&,points,{-2}],{-3}],Map[If[vector[[2]] ≥ 0,ToSequence[],Polygon[#]]&,Map[#+vector*(0-#[[2]])/vector[[2]]&,points,{-2}],{-3}]}]], {GrayLevel[0],{Line[{{0,0,0},{x,0,0},{x,y,0},{0,y,0},{0,0,0}}],Line[{{0,0,z},{x,0,z},{x,y,z},{0,y,z},{0,0,z}}],Line[{{0,0,0},{0,0,z}}],Line[{{x,0,0},{x,0,z}}],Line[{{x,y,0},{x,y,z}}],Line[{{0,y,0},{0,y,z}}]}}, ToSequence @@ epilog},ToSequence @@DeleteCases[{opts},(Beam->_) | (ClipDistance->_) | (Direction->_) | (Epilog->_) | (Orientation->_) | (Prolog->_) | (RoomDimensions->_) | (RoomLength->_)| (RoomWidth->_) | (WindowPoints->_),{1}],Lighting->{{"Ambient",White}}, Boxed->False,PlotRange->{{-clip,x+clip},{-clip,y+clip},{-clip,z+clip}}]
]

Explanation

This notebook requires internet connectivity to access astronomical data.

Sunlight[t,opts] creates a graphic which illustrates the regions of direct sunlight in a rectangular room at time t.  The time t should be expressed in the same format as the output of the Mathematica command DateList.  The geometry of the room is determined by the options opts.  A Cartesian coordinate system is superimposed on the room, with the floor lying in the x-y plane.  One wall is assumed to lie in the x-z plane, and another wall is assumed to lie in the y-z plane.  The length l of the room along the x-axis and width w of the room along the y-axis are specified with the option  RoomDimensions->{l,w}.  The shape and placement of a window is determined by WindowPoints->pts, where pts is a list of Cartesian coordinates which define the perimeter of the window.  The orientation of the room is determined by Orientation->θ  Degree, where θ is the number of degrees away from north (measured in a clockwise direction) that one is facing if one stands parallel to the x-axis and faces in the positive direction.  The latitude and longitude (and, if desired, the elevation) of the room should be specified by setting the Mathematica $GeoLocation variable.

$GeoLocation={40.11,-88.18};

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}}]

These points correspond to a semicircular window.

semicircle =Table[{0,9+x,7+Sqrt[4^2-(x-4)^2]},{x,0,8}]

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->semicircle]

Multiple windows may also be specified.  If list1, list2, etc. are the lists of points which define multiple windows, then use WindowPoints->{list1,list2,…} to represent those windows.

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle}]

The windows need not be on the same wall.

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle,{{5.5,26,1}, {5.5,26,8},{13.5,26,8},{13.5,26,1}}}]

The following is an animation of a room on 21 December 2000.  Each frame represents 6 minutes.

Manipulate[Sunlight[DateList[{2000,12,21,hour},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle,{{5.5,26,1}, {5.5,26,8},{13.5,26,8},{13.5,26,1}}}],{hour,7.25,16.25,.1}]

The beams of sunlight can be displayed with the option Beam->True.

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle,{{5.5,26,1}, {5.5,26,8},{13.5,26,8},{13.5,26,1}}}, Beam->True]

Mathematica graphics primitives prms can be included with the options Prolog->{prms} and Epilog->{prms}.  Primitives specified by Prolog are rendered before all other primitives, and primitives specified by Epilog are rendered after all other primitives.  Prolog and Epilog can be used to simulate objects in the room.  Note that these simulated objects cannot cast shadows or interact with the beams of light, though.

The image below includes a simulated desk.

Sunlight[DateList[{2000,12,21,8},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle,{{5.5,26,1}, {5.5,26,8},{13.5,26,8},{13.5,26,1}}},Epilog->{Cuboid[{0,4,0},{8,7,2.5}]}]

An animation can be made which illustrates how the beams of sunlight intersect the desk throughout the course of the day on 21 December 2000.

Manipulate[Sunlight[DateList[{2000,12,21,hour},TimeZone->-6], Orientation->275 Degree,RoomDimensions->{19,26},WindowPoints->{{{0,9,1},{0,9,6.5},{0,17,6.5},{0,17,1}},semicircle,{{5.5,26,1}, {5.5,26,8},{13.5,26,8},{13.5,26,1}}},Epilog->{Cuboid[{0,4,0},{8,7,2.5}]},Beam->True],{hour,7.25,16.25,.1}]