Simulating Sunlight

December 1999

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

Definitions

<< Astronomer`HomeSite`;

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, DaylightSavingTime -> Automatic, 
          Orientation -> 0 Degree, Epilog -> {}, 
          Latitude -> 40.1 Degree, Longitude -> -88.2 Degree, 
          Prolog -> {}, RoomDimensions -> {1, 1}, 
          TimeZone -> -6, 
          WindowPoints -> {{0, 0, 0}, {0, 0, 1}, {0, 1, 1}, {0, 1, 
                0}}};

Sunlight[{month_, day_, hour_}, 
        opts___]  := 
Module[{beam, clip, daylight, dir, 
          epilog, floor, lat, len, long, orient, points, prolog, vector, wid, 
          x, y, z, 
          zone}, 
beam  =  Beam  /.  {opts}  /.  
            Options[Sunlight]; 
clip  =  ClipDistance  /.  {opts}  /. 
            Options[Sunlight]; 
daylight  =  DaylightSavingTime  /.  {opts}  /. 
            Options[Sunlight]; 
epilog  =  Epilog  /.  {opts}  /. 
            Options[Sunlight]; 
lat  =  Latitude  /.  {opts}  /. 
            Options[Sunlight]; 
long  =  Longitude  /.  {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]; 
zone  =  TimeZone  /.  {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]; 
SetLocation[GeoLatitude -> lat, 
          GeoLongitude -> long, GeoAltitude -> 0.0, 
          TimeZone -> zone]; 
z  = 
          Max[Flatten[
              Map[#[[3]] &,  points, {-2}]]]; 
If[
          daylight === Automatic, 
          If[ModifiedJulianDay[{2000, 4, 2, 3, 0, 0}] ≤ 
              ModifiedJulianDay[{2000, month, day, hour, 0, 
                  0}]  ≤ 
              ModifiedJulianDay[{2000, 10, 29, 2, 0, 0}], daylight = True, 
            daylight = False]]; 
vector = 
          N[-{Cos[orient - Azimuth]*Cos[Altitude], 
                  Sin[orient  - Azimuth]*Cos[Altitude], 
                  Sin[Altitude]} /.   HorizonCoordinates[
                Sun, {2000, month, day, hour - If[daylight === True, 1, 0], 
                  0, 0}]]; 
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 -> _)  |  (DaylightSavingTime -> _)  |  (Direction -> _)  |  (Epilog -> _)  |  (Latitude -> _)  |  (Longitude -> _)  |  (Orientation -> _)  |  (Prolog -> _)  |  (RoomDimensions -> _)  |  (RoomLength -> _) |  (RoomWidth -> _) |  (TimeZone -> _)  |  (WindowPoints -> _), {1}], AmbientLight -> RGBColor[1, 1, 1], 
          LightSources -> {},  Boxed -> False, 
          PlotRange -> {{-clip, x + clip}, {-clip, 
                y + clip}, {-clip, 
                z + clip}}]
]

Explanation

This notebook requires the Scientific Astronomer package.

Sunlight[{m,d,h},opts] creates a graphics object which illustrates the regions of direct sunlight in a rectangular room on month m, day d, and hour h of the year 2000.  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 lat, longitude lng, and time zone zn may be specified with the options Latitude->lat, Longitude->lng, and TimeZone->zn.  If these options are omitted, then the latitude and longitude are assumed to be that of Champaign, IL, and the time zone is assumed to be -6, which is also standard for Champaign.  The default option DaylightSavingTime->Automatic specifies that Daylight Saving Time (as practiced in most of the United States) is automatically taken into account.  The option DaylightSavingTime->False disables Daylight Saving Time and DaylightSavingTime->True forces Daylight Saving Time

Show[
        Sunlight[{12, 21, 8},  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}]

Show[
        Sunlight[{12, 21, 8},  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.

Show[
        Sunlight[{12, 21, 8},  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.

Show[
        Sunlight[{12, 21, 8},  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.  Each frame represents 6 minutes.

Do[
        Show[Sunlight[{12, 21, i},  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}}}]], {i, 7.25, 
          16.25,  .1}];

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

Show[
        Sunlight[{12, 21, 8},  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.

Show[
        Sunlight[{12, 21, 8},  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.

Do[
        Show[Sunlight[{12, 21, i},  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]], {i, 7.25, 16.25,  .1}];