The Labs.Com Issue_06_Sundials

TPJ: Issue_06_Sundials

This is a collection of programs published by The Perl Journal. You can download all source-code also from TPJ: Programs.

1. sundial.pl
2. sd.dat
3. sd.ps
4. More Samples on Sundials

Download sundial.pl

#!/usr/bin/perl

#

# A Sundial Generator

#

# (c) J.L. Redford

#

$VersionDate = "03/31/97";

#------------------- Main code -----------------------

&InitUserParameters;

&InitGlobals;

&StartPage ($OutputFileName);

&Draw2HorizBorder;

&DrawDial (1);

&DrawDial (0);

&ClosePage;

exit;

#-------------- User Interface ----------------------

sub InitUserArrays {

# Table of time zone names versus the center longitude of the zone

%CentralLongitude = (

"EST", 75,

"CST", 90,

"MST", 105,

"PST", 120

);

$AllowedTimeZones = join ('|', keys (%CentralLongitude));

$AllowedMottos = join ('|', (

"I Mark Time, but Time Marketh Thee",

"Tempus Fugit",

"I Count Only Sunny Hours",

"[your motto here]"));

# Defines things that the user can change.  The key

# field is the variable name.  In the other field are the

# default value, a description string, a code indicating the

# parameter type, and a range of allowed values

%UserParameters = (

"OutputFileName", "1|sd.ps|Postscript output file name|STR",

"Latitude", "2|42:20|Latitude of dial in degrees north of the equator|\

NUM|-90|90",

"Longitude", "3|71:00|Longitude of dial in degrees west of Greenwich|N\

UM|-180|180",

"TimeZone", "4|EST|Time Zone of dial|STR|$AllowedTimeZones",

"LastWinterTime","5|15:10|Latest time to show on the winter dial|HM|15\

:00|18:00",

"FirstSummerTime","6|5:10|Earliest time to show on the summer dial|HM|\

3:00|9:00",

"Motto", "7|Tempus Fugit|Sentitious saying for decoration|STR|",

);

}

sub Prompt {

local ($var) = @_;

($pos,$default,$desc,$type,@Range) = split ('\|',$UserParameters{$va\

r});

print "$var [$default] ";

$response = <>;

chop $response;

if ($response =~ /\?|help/) {

print "$desc\n";

print "$var [$default] ";

$response = <>;

chop $response;

if ($response =~ /\?|help/) {

print "using default value\n";

$response = $default;

}

}

elsif ($response eq "") {

$response = $default;

}

eval ("\$$var = \'$response\'");

if ($var eq "OutputFileName") {

($root) = $response =~ /(\w+)\.?\w*$/;

$DatFileName = "$root.dat";

open (DAT, ">".$DatFileName) ||

die "Cannot open file to record parameters, $DatFileName";

open (PS, ">".$OutputFileName) ||

die "Cannot open Postscript output file, $OutputFileName";

print PS "%\n";

print PS "% Parameters used for this dial:\n";

}

print DAT "  # $desc\n";

print DAT "\$$var = \"$response\";\n";

print PS "% $var = $response\n";

}

sub InitUserParameters {

&InitUserArrays;

if ($ARGV[0] =~ /\?|help/) {

print <<End_of_usage;

End_of_usage

exit;

}

elsif ($ARGV[0] ne "") {

require $ARGV[0];

}

else {

print "Sundial Generator $VersionDate (c) J.L. Redford\n";

print "[]s show default values.\n";

print "Type ? or help to get a parameter description\n";

foreach $v (keys %UserParameters) {

($pos,@foo) = split ('\|', $UserParameters{$v});

$Params {$pos} = $v;

}

foreach $v (sort keys %Params) {

&Prompt ($Params {$v});

}

}

}

#-------------- Sundial Mathematics -----------------

# Return a 3-vector of the sun's position

sub GetSunV {

local ($d, $t) = @_;

# Find the day of the year as an integer

$day = int ($d * $days_per_year / $twopi) + $VernalEquinox;

if ($day < 0) { $day += 365; }

# Look up the equation of time for that day and add it to the

# clock time.

$t += ($EQT[$day]) * $twopi / (60*60*24);

# Correct for longitude

$t += $LongitudeCorrection;

#  return ( &RotX ($l, &RotZ (-$t-$d, &RotX ($e, &RotZ ($d, (1,0,0))))\

))

# Rotate by dial strike and dip, by latitude, by time, and by

# sun's position in the year

return ( &RotX ($l, &RotZ (-$t-$twopi/4, &RotX ($e * sin ($d), (0,1,\

0)))));

}

# Return a 2-vector of the tip of the gnomon shadow for a gnomon at (0\

,0,1)

sub GetTipV {

local (@S);

@S = &GetSunV ($_[0], $_[1]);

$_[0] * 360/$twopi, $_[1] * 360/$twopi, $S[0], $S[1], $S[2];

return (($S[2] > 0) ? (- $S[0]/$S[2], - $S[1]/$S[2]) : 0);

}

# Rotate a vector @in by $a radians around the Z axis

sub RotZ {

($a,@in) = @_;

return (

(cos ($a) * $in[0] - sin ($a) * $in[1]),

(sin ($a) * $in[0] + cos ($a) * $in[1]),

$in[2]

);

}

# Rotate a vector @in by $a radians around the X axis

sub RotX {

($a,@in) = @_;

return (

$in[0],

(cos ($a) * $in[1] - sin ($a) * $in[2]),

(sin ($a) * $in[1] + cos ($a) * $in[2])

);

}

sub InitGlobals {

$twopi = 2 * 3.1415926536;

$days_per_year = 365.24;

# day on which the Sun is on the celestial equator

$VernalEquinox = 78;

# tilt of the earth's axis with respect to the celestial equator

$e = 23.5 * $twopi / 360;

# Convert latitude in degrees:minutes to a decimal

($deg,$min) = split (':',$Latitude);

$LatitudeDecimal = &HourMinToDecimal ($Latitude);

# Amount to rotate dial by to account for the dial latitude

$l = (90 - $LatitudeDecimal) * $twopi / 360;

$LongitudeDecimal = &HourMinToDecimal ($Longitude);

$LongitudeCorrection =

($CentralLongitude {$TimeZone} - $LongitudeDecimal) * $twopi / 360\

;

# Latitude and Longitude inscription

$LatDeg = int $LatitudeDecimal;

$LatMin = int (($LatitudeDecimal - $LatDeg)*60);

if ($LatDeg > 0) {

$LatLongInscription = sprintf ("%d\\312%2d\'N ", $LatDeg, $LatMin)\

;

}

else {

$LatDeg = - $LatDeg; $LatMin = -$LatMin;

$LatLongInscription = sprintf ("%d\\312%2d\'S ", $LatDeg, $LatMin)\

;

}

$LongDeg = int $LongitudeDecimal;

$LongMin = int (($LongitudeDecimal - $LongDeg)*60);

if ($LongDeg >= 0 && $LongDeg <= 180) {

$LatLongInscription .= sprintf ("%d\\312%2d\'W", $LongDeg, $LongMi\

n);

}

elsif ($LongDeg > 180) {

$LongDeg = 360 - $LongDeg; $LongMin = 60 - $LongMin;

$LatLongInscription .= sprintf ("%d\\312%2d\'E", $LongDeg, $LongMi\

n);

}

elsif ($LongDeg < 0) {

$LongDeg = - $LongDeg; $LongMin = -$LongMin;

$LatLongInscription .= sprintf ("%d\\312%2d\'E", $LongDeg, $LongMi\

n);

}

$LastWinterTimeDec = &HourMinToDecimal ($LastWinterTime);

$FirstSummerTimeDec = &HourMinToDecimal ($FirstSummerTime);

# Translate hour names to Roman Numerals.  Notice that the hours

# past noon wrap around.

%RomanNumerals = (

1, "I",    2, "II",    3, "III",    4, "IV",

5, "V",    6, "VI",    7, "VII",    8, "VIII",

9, "IX",  10, "X",    11, "XI",    12, "XII",

13, "I",  14, "II",   15, "III",   16, "IV",

17, "V",  18, "VI",   19, "VII",   20, "VIII",

21, "IX", 22, "X",    23, "XI",    24, "XII"  );

# Equation of time in seconds for each day of the year, starting Jan\

1

@EQT = (

-188,-216,-245,-274,-303,-327,-352,

-377,-401,-426,-451,-475,-500,-525,

-550,-568,-586,-604,-622,-641,-659,

-677,-695,-713,-732,-741,-751,-761,

-771,-781,-791,-801,-811,-821,-831,

-841,-842,-844,-845,-847,-848,-850,

-851,-853,-854,-856,-850,-844,-838,

-832,-827,-821,-815,-809,-803,-797,

-786,-774,-763,-751,-739,-728,-716,

-704,-689,-674,-659,-644,-629,-613,

-598,-583,-568,-553,-535,-517,-499,

-482,-464,-446,-429,-411,-393,-376,

-357,-339,-321,-303,-285,-267,-249,

-231,-213,-195,-176,-157,-138,-119,

-100, -81, -62, -43, -24,  -5,  14,

24,  34,  44,  54,  65,  75,  85,

95, 105, 116, 124, 132, 140, 148,

157, 165, 173, 181, 189, 198, 200,

203, 205, 208, 211, 213, 216, 218,

221, 224, 221, 218, 215, 212, 210,

207, 204, 201, 198, 196, 187, 179,

171, 163, 155, 146, 138, 130, 122,

114, 106,  94,  82,  71,  59,  48,

36,  24,  13,   1, -10, -23, -36,

-49, -62, -75, -88,-101,-114,-127,

-140,-151,-163,-175,-187,-199,-211,

-223,-235,-247,-259,-267,-276,-285,

-293,-302,-311,-319,-328,-337,-346,

-349,-353,-357,-361,-365,-368,-372,

-376,-380,-384,-381,-379,-377,-374,

-372,-370,-368,-365,-363,-361,-359,

-350,-341,-333,-324,-316,-307,-298,

-290,-281,-273,-259,-245,-231,-217,

-203,-189,-175,-161,-147,-134,-115,

-97, -79, -61, -43, -25,  -7,  10,

28,  46,  65,  85, 106, 127, 147,

168, 189, 209, 230, 251, 272, 293,

314, 335, 356, 377, 399, 420, 441,

462, 483, 503, 523, 542, 562, 581,

601, 621, 640, 660, 679, 696, 712,

728, 744, 760, 776, 792, 808, 824,

841, 851, 862, 872, 883, 894, 904,

915, 925, 936, 947, 950, 953, 956,

959, 962, 966, 969, 972, 975, 978,

981, 976, 971, 965, 960, 954, 949,

944, 938, 933, 927, 914, 900, 886,

873, 859, 845, 832, 818, 804, 790,

769, 748, 727, 705, 684, 663, 641,

620, 599, 578, 551, 524, 497, 470,

443, 416, 389, 362, 335, 308, 279,

249, 220, 190, 161, 131, 101,  72,

42,  13, -15, -44, -73,-101,-130,

-159,

);

}

# Convert a string of the form hours:minutes to a decimal

sub HourMinToDecimal {

local ($hour,$min) = split (':', $_[0]);

return ($hour + $min / 60);

}

# Convert a number to a string of the form hours:minutes

sub DecimalToHourMin {

local ($hour,$min);

$hour = int $_[0];

$min = int (($_[0] - $hour) * 60);

return ("$hour:$min");

}

#----------------- Drawing subroutines -----------------------

# Draw the non-dial parts of the two horizontal dials

sub Draw2HorizBorder {

print PS "-4 -5.25 4 5.25 set_corners\n";

print PS "18 selectsize\n";

print PS "heavy setlinewidth\n";

&DrawCenteredString (0,0,$Motto);

&DrawCenteredString (0,-4.75, $LatLongInscription);

#  &StartLine (-4, -5.25);

#  &LineTo (-4, 5.25);

#  &LineTo (4, 5.25);

#  &LineTo (4, -5.25);

#  &LineTo (-4, -5.25); &DrawLine;

}

# Draw a spring or winter dial

sub DrawDial {

local ($SpringDial) = @_;

# Find the bounding box so that the point for

# LastWinterTime is in the lower left and FirstSummerTime is

# in the upper right.

@T1 = &GetTipV (-$twopi/4, &HoursToRadians ($LastWinterTimeDec));

@T2 = &GetTipV ($twopi/4, &HoursToRadians ($FirstSummerTimeDec));

if ($SpringDial) {

&FillHalf (1, @T1, -$T1[0], $T2[1]);

}

else {

&FillHalf (0, -$T2[0], $T1[1], @T2);

}

$StartHour = 5;

$EndHour = 19;

$StartDate = $SpringDial ? -$twopi/4 : $twopi/4;

$EndDate = $SpringDial ? $twopi/4 : $twopi*3/4;

$Legend = $SpringDial ? "Winter and Spring" : "Summer and Fall";

print PS "heavy setlinewidth\n";

print PS "0 0 mark\n";

print PS "18 selectsize\n";

&DrawCenteredString (0,1,$Legend);

&StartLine (-0.5, 0.5); &LineTo (0.5, 0.5); &DrawLine;

print PS "14 selectsize\n";

for ($hour = $StartHour; $hour <= $EndHour; $hour += 1) {

$t = &HoursToRadians ($hour);

&LabelHour ($hour);

$DashedLine = 0;

print PS "heavy setlinewidth\n";

&DrawConstantLine ($StartDate, $EndDate, $twopi/200, 1);

$DashedLine = 1;

print PS "light setlinewidth\n";

for ($min = 10; $min < 60; $min += 10) {

$t = &HoursToRadians ($hour + $min/60);

&DrawConstantLine ($StartDate, $EndDate, $twopi/200, 1);

}

}

$DashedLine = 0;

print PS "heavy setlinewidth\n";

foreach $d ($StartDate, $EndDate) {

&DrawConstantLine (

&HoursToRadians($StartHour),

&HoursToRadians($EndHour),

(10/60) * $twopi / 24, 0);

}

#  &CloseImage;

}

# Label a constant-time line with its time in Roman numerals

sub LabelHour {

local ($hour) = @_;

local ($t1, $t2, $d, @LP1, @LP2, $dx, $slope);

$Inside = ($hour > 15) || ($hour < 9);

$t1 = &HoursToRadians ($hour);

$t2 = &HoursToRadians ($hour+0.01);

$d = $Inside ? ($twopi/4) : (-$twopi/4);

@LP1 = &GetTipV ($d, $t1);

@LP2 = &GetTipV ($d, $t2);

$dx = $LP1[0] - $LP2[0];

if ($dx == 0) {$dx = 0.000001;}

$slope = ($LP1[1] - $LP2[1]) / $dx;

if ($slope > 1000) { $slope = 1000; }

if ($slope < -1000) { $slope = -1000; }

$function = $Inside ? "ce" : "cet";

printf PS ("%.3f %.3f %.3f (%s) %s\n",  $slope, $LP1[0], $LP1[1],

$RomanNumerals {$hour}, $function);

}

# Convert hours (in 24 hour time) to radians

sub HoursToRadians {

return (($_[0] - 18) * $twopi / 24);

}

# Draw a line of constant date or constant time

sub DrawConstantLine {

local ($Start, $End, $Incr, $ConstantTime) = @_;

$first = 1; $inline = 0; $inp = 1;

for ($i = $Start; $i <= $End; $i += $Incr) {

@T = $ConstantTime ? &GetTipV ($i, $t) : &GetTipV ($d, $i);

if ($#T) {

$in = &InsideFrame (@T);

if ($in) {

if (! $inp) { &StartLine (@Tp); &LineTo (@T); $inline+=2;}

elsif ($first) { &StartLine (@T); $inline++;}

elsif ($inline) { &LineTo (@T); $inline ++;}

$first = 0;

}

else {

if ($inp) { &LineTo (@T); $inline++; }

}

@Tp = @T; $inp = $in;

}

}

if ($inline > 1) { &DrawLine; }

}

# Return 1 if the input vector is within the frame defined by

# LeftEdge, RightEdge, TopEdge, and BottomEdge

sub InsideFrame {

return (($_[0] < $RightEdge && $_[0] > $LeftEdge) &&

($_[1] < $TopEdge   && $_[1] > $BottomEdge));

}

#----- Postscript Graphics Functions -----------------

sub StartLine { $PSLine = sprintf ("%.3f %.3f m ",$_[0], $_[1]); }

sub LineTo {    $PSLine .= sprintf ("%.3f %.3f l ", $_[0], $_[1]); }

sub DrawLine {

if ($DashedLine) {

# Change the odd-numbered "l"s to "f"s, and the evens to "m"s

$PSLine =~ s/(\S+ \S+) l (\S+ \S+) l/\1 f \2 m/g;

print PS $PSLine,"\n";

}

else {

$PSLine =~ s/l $/f\n/;

print PS $PSLine;

}

}

sub DrawCenteredString {

printf PS ("%.3f %.3f (%s) cshow\n", @_);

}

sub FillHalf {

$bool = $_[0] ? "true" : "false";

printf PS ("%s %.3f %.3f %.3f %.3f place\n", ($bool, @_));

$LeftEdge = $_[1];

$BottomEdge = $_[2];

$RightEdge = $_[3];

$TopEdge = $_[4];

}

sub StartPage {

print PS <<EOF;

/m {moveto} def

/l {lineto} def

/s {stroke} def

/f {lineto stroke} def  % finish a line

/makestr {/sname xdef dup 20 string cvs sname xdef} def

/heavy 0.02 def         % for heavy lines

/light 0.01 def         % for light lines

/xdef {exch def} def        % saves text

% select a font size given the size in points

/selectsize {       % usage: size selectsize

/size xdef

/Times-Roman findfont

1 dup idtransform pop     % get size of one point

75 mul 24 div            % oddball scaling

dup /point_size xdef        % store for use in margins

size mul scalefont setfont

gsave

newpath

0 0 moveto

(A) true charpath                % Use an "A" for the height of al\

l strings.

flattenpath pathbbox

4 1 roll pop pop pop

grestore

/sh xdef

} def

/set_corners {              %usage: llx lly urx ury set_corners

/new_ury xdef       % pull the args off the stack and name them

/new_urx xdef

/new_lly xdef

/new_llx xdef

clippath pathbbox       % get coords of present bounding box

/ury xdef

/urx xdef

/lly xdef

/llx xdef

% scale so that the y direction goes from top edge to bottom edge

/scale_factor ury lly sub new_ury new_lly sub div def

% get x translation

new_urx scale_factor mul llx add

% then y translation

ury lly sub 2 div lly add

translate

scale_factor dup scale

} def

% Place a string so that its bottom center is as close as possible to \

(x,y) but

% doesn't cross the line defined by slope and (x,y)

/ce                 % usage: slope x y (s) ce

{

/s xdef

s stringwidth pop % get string width

point_size 4 mul add /sw xdef       % including a 2-point margin

point_size 4 mul sh add /shm xdef

/cx sw 2 div def

/cy shm 2 div neg def

center

} def

% Place a string so that its top center is as close as possible to (x,\

y) but

% doesn't cross the line defined by slope and (x,y)

/cet                 % usage: slope x y (s) cet

{

/s xdef

s stringwidth pop  % get string width

point_size 4 mul add /sw xdef       % including a 2-point margin

point_size 4 mul sh add /shm xdef

/cx sw 2 div neg def

/cy shm 2 div def

center

} def

/center {

/y xdef

/x xdef

/slope xdef

slope 0 lt {/cx cx neg def} if      % negate cx if slope negative

/dx slope cy mul slope dup mul cx mul sub slope dup mul 1 add div de\

f

/dy slope cx dx add mul cy sub def

x dx add sw 2 div sub point_size 2 mul add

y dy add sh 2 div sub

moveto

s show

} def

% Draw a centered string

/cshow {

/s xdef

/y xdef

/x xdef

x s stringwidth pop 2 div sub y moveto s show

} def

% set the transform so as to fill the upper or lower part of a sheet

/place {              %usage: upper_or_lower llx lly urx ury place

/new_ury xdef       % pull the args off the stack and name them

/new_urx xdef

/new_lly xdef

/new_llx xdef

/upper   xdef

clippath pathbbox       % get coords of present bounding box

/ury xdef

/urx xdef

/lly xdef

/llx xdef

% scale so that the x direction goes from left edge to right edge

/scale_factor urx llx sub new_urx new_llx sub div def

% get x translation

new_llx neg scale_factor mul llx add

% then y translation

ury lly sub 2 div lly add       % find midline

upper

{new_lly neg scale_factor mul add}

{new_ury scale_factor mul sub}

ifelse

translate

scale_factor dup scale

% orient so that long side is along paper

} def

% draw a little marker at (x,y)

/mark {      % usage: x y mark

/y xdef

/x xdef

x 0.05 add y moveto x 0.05 sub y lineto stroke  % draw horizontal

x y 0.05 add moveto x y 0.05 sub lineto stroke  % draw vertical

x y 0.05 0 360 arc stroke                       % draw circle

} def

/d1_font_size 1 def

EOF

}

sub ClosePage {

print PS "showpage\n";

close PS;

}

Issue_06_Sundials 2. sd.dat

sd.dat

Issue_06_Sundials 3. sd.ps

sd.ps

Issue_06_Sundials 4. More Samples on Sundials

• Issue_06_Sundials