#!/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|NUM|-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{$var});
  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, $LongMin);
  }
  elsif ($LongDeg > 180) {
    $LongDeg = 360 - $LongDeg; $LongMin = 60 - $LongMin;
    $LatLongInscription .= sprintf ("%d\\312%2d\'E", $LongDeg, $LongMin);
  }
  elsif ($LongDeg < 0) {
    $LongDeg = - $LongDeg; $LongMin = -$LongMin;
    $LatLongInscription .= sprintf ("%d\\312%2d\'E", $LongDeg, $LongMin);
  }

  $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 all 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 def
  /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;
}