function [x,y] = sunshadow(latitude, day, hour)
%SUNSHADOW Plot the shadow of a unit-length pin for a given time and place.
% This function is useful for creating sundials.
% [x,y] = sunshadow(latitude) returns the location of the tip of the shadow
% of a sundial gnomon (pin). The pin is assumed to be unit length and
% perpendicular to the plane of the sundial's surface.
% [x,y] = sunshadow(latitude, day) returns the location of the shadow tip
% for a given day number (or list of day numbers)
% [x,y] = sunshadow(latitude, day, hour) returns the location of the shadow
% tip for a give day number and hour number.
% NOTE: Only the day number can be vectorized
%
% Example
% sunshadow(42)
%
% Reference: The algorithms in this file were adapted from
% http://www.de-zonnewijzerkring.nl/eng/
% Copyright 2009 The MathWorks, Inc.
% Definitions
% phi
% latitude of the place of the dial. -90 <= phi <= 90,
% northern latitude positive, southern latitude negative.
% g
% length of (pin) gnomon perpendicular to the dial's plane.
% The tip of the gnomon is the shadow casting point. g > 0.
% i
% inclination of the plane: the zenith distance of the
% endpoint of the (pin) gnomon. 0 <= i < 180.
% d
% declination of the plane: azimuth of the gnomon: south = 0,
% positive to west, negative to east. -180 <= d <= 180.
% E
% equation of time in degrees. In November E is positive.
% Conversion minutes of time into degrees: 1 minute = 0.25 degrees.
% decl
% declination of the sun, pos. if sun in northern hemisphere,
% neg. in southern hemisphere. -23.5 <= decl <= 23.5
% t
% hourangle of the sun: noon = 0, positive. to west, negative
% to east. -180 <= t <= 180.
% x,y
% coordinates of shadow point. Footpoint gnomon = 0,0. x to right,
% y upwards. For horizontal dial x east, y north.
% dn
% daynumber. January first at 00:00:00 is daynumber 1.0 , at
% 12:00:00 is daynumber 1.5.
% ...
% diverse other variables are used in the procedures as x0,y0,z0,
% x1, y1, z1, day, mth and so on. See the procedures.
% Unused in this file
% SM
% standard meridian of time zone. Greenwich = 0, positive
% to west, negative to east. -180 <= SM <= 180.
% LM
% local meridian. (longitude of the place of the sundial.)
% Greenwich = 0, positive to west, negative to east. -180 <= LM <= 180.
% LC
% longitude correction LC = SM - LM.
% v
% height of a style, parallel to the earth axis. -90 <= v <= 90.
% b
% angle of substyle, measured from y-axis. pos. anti clockwise.
% -180 <= b <= 180.
% ts
% hourangle of substyle. -180 <= ts <= 180.
%
% While these definitions are, of course, arbitrary, many of them are
% common in gnomonics and the others are constrained by the requirement
% that the method should be universal all over the world. If you want to
% change the definitions , you may also have to change some of the
% formulae, or introduce a conversion sub-routine.
%
% English-speaking people nowadays use the word gnomon to refer to the
% whole of the line casting the shadow, but in this paper gnomon is used
% in its older and more traditional sense. Here the gnomon is a pin,
% perpendicular to the plane.
% Default location is Boston, MA
% North Latitude 42.36 deg West Longitude 71.19 deg
if nargin < 1
latitude = 42.36;
end
if nargin < 2
% Use today's day number
dv = datevec(now);
day = floor(now - datenum(dv(1),1,1)) + 1;
end
if nargin < 3
hour = 12;
end
% d2r = degrees to radians conversion
d2r = pi/180;
phi = latitude*d2r;
% The gnomon is of length 1, pointing straight up on a flat surface
g = 1;
i = 0*d2r;
d = 0*d2r;
dn = day;
u = hour;
% Calculate decl for daynumber ( dn + 0.5 )
% Calculate equation of time E for daynumber ( dn + 0.5 )
[decl, E] = eqnoftime(dn+0.5);
% Calculate hourangle
% E is in seconds of time and u is in hours
% Two conversion factors are applied to make t come out in radians
% E is multiplied by (1 min/60 sec)*(1 hr/60 min)*(360 deg/24 hr)
t = ((u - 12)*15 + E*15/(60*60))*d2r;
% Call mainprocedure
[x, y] = mainproc(phi, g, i, d, decl, t);
if nargout == 0
plot([0 x],[0 y])
axis equal
end
end
% =========================================================================
function [x, y] = mainproc(phi, g, i, d, decl, t)
% The Main Procedure
% The aim of the main procedure is to convert a given sun's position,
% defined by the angle of declination decl and the sun's
% hourangle t into the coordinates x, y of the shadow point of the (pin)
% gnomon or to output 'point isn't real'.
% Just simply calculate the values x0,y0,z0, x1,y1,z1, x2,y2,z2,
% x3,y3,z3, and finaly x,y as in the 5 routines below.
% In these routines you see 2 decision points to find out if a point is or isn't real.
% In : phi, g, i, d, decl, t
% Out : x, y or 'point isn't real'
x0 = sin(t) .* cos(decl);
% translation decl, t into x0, y0, z0
y0 = cos(t) .* cos(decl);
z0 = sin(decl);
R = pi/2 - phi;
% translation x0, y0, z0, into x1, y1, z1
x1 = x0;
% by rotation around x-axis through angle 90 - phi
y1 = y0 * cos(R) - z0 * sin(R);
z1 = y0 * sin(R) + z0 * cos(R);
% if z1 < 0 point isn't real: sun is beneath the horizon.
if z1<0
x=NaN;
y=NaN;
return
end
R = d;
% translation x1, y1, z1, into x2, y2, z2
x2 = x1 * cos(R) - y1 * sin(R);
% by rotation around z-axis through angle d
y2 = x1 * sin(R) + y1 * cos(R);
z2 = z1;
R = i;
% translation x2, y2, z2, into x3, y3, z3
x3 = x2;
% by rotation around x-axis through angle i
y3 = y2 * cos(R) - z2 * sin(R);
z3 = y2 * sin(R) + z2 * cos(R);
% if z3 <= 0 point isn't real: sun isn't above the dial.
if z3 <= 0
x = NaN;
y = NaN;
return
end
x = x3 .* (g./z3);
% translation x3, y3, z3, into x, y
y = y3 .* (g./z3);
% these are the wanted coordinates of the shadow point.
end
% =========================================================================
function [decl, E] = eqnoftime(dayNr)
% Equation of time calculations
% d2r = degrees to radians conversion
d2r = pi/180;
L = (dayNr*360/365.2422 - 80.535132) * d2r; % radians
% Equation of time
E = - 107.0605*sin( L) - 428.6697*cos( L) ...
+ 596.1009*sin(2*L) - 2.0898*cos(2*L) ...
+ 4.4173*sin(3*L) + 19.2776*cos(3*L) ...
- 12.7338*sin(4*L); % seconds (of time)
lambda = L + (0.4277*sin( L) + 1.8664*cos( L) ...
+ 0.0180*sin(2*L) + 0.0087*cos(2*L)) * d2r; % radians
epsilon = 23.43954 * d2r; % radians
decl = asin(sin(lambda) * sin(epsilon)); % radians
end
