pom

#!/usr/bin/perl

=begin metadata

Name: pom
Description: display the phase of the moon
Author: Rocco Caputo, troc@netrus.net
License: perl

=end metadata

=cut


use strict;

use File::Basename qw(basename);
use Getopt::Std qw(getopts);
use Math::Trig qw(atan deg2rad rad2deg tan);
use POSIX qw(floor fmod);

use constant EX_SUCCESS => 0;
use constant EX_FAILURE => 1;

use constant LINES_DEFAULT => 25;
use constant COLS_DEFAULT  => 80;
use constant ASPECT_RATIO => 5 / 7;

my $Program = basename($0);

{ my $vt100_compatible = 0;
  $vt100_compatible ||= $ENV{TERM}    =~ /vt100|xterm|ansi/i;
  $vt100_compatible ||= $ENV{TERMCAP} =~ /vt100|xterm|ansi/i;
  if ($^O eq 'VMS') {
    $vt100_compatible = grep(/VT\d{2,}/i, `SHOW TERMINAL $ENV{TT}`)
  }
  my $cursor_home  = $vt100_compatible ? "\e[0;0H" : ('-' x 78 . "\n");
  my $screen_clear = $vt100_compatible ? "\e[2J" : '';
  sub CURSOR_HOME  { $cursor_home }
  sub SCREEN_CLEAR { $screen_clear }
}

#
## Astronomical constants

sub EPOCH    () { 2444238.5    } # 1980 January 0.0
sub ELONGE   () { 278.833540   } # ecliptic longitude of the Sun at EPOCH
sub ELONGP   () { 282.596403   } # ecliptic longitude of the Sun at perigee
sub ECCENT   () {   0.01671542 } # Earth's orbit's eccentricity
sub MMLONG   () {  64.975464   } # moon's mean longitude at EPOCH
sub MMLONGP  () { 349.383063   } # mean longitude of the perigee at EPOCH

#
## Helper functions

sub fixangle { fmod(($_[0] - (360 * (floor($_[0] / 360)))), 360); }

sub usage {
  warn "usage: $Program [-d] [-e] [[[[[[cc]yy]mm]dd]HH]]\n";
  exit EX_FAILURE;
}

sub checknum {
  my $n = shift;
  my $lab = shift || 'number';
  if ($n !~ m/\A[0-9]+\Z/) {
    warn "$Program: bad $lab: '$n'\n";
    exit EX_FAILURE;
  }
  return int($n);
}

#
## Parse the command line, filling in with the current GMT

my ($cc, $yy, $mm, $dd, $hh);
my ($gm_yy, $gm_mm, $gm_dd, $gm_hh) = (gmtime(time))[5,4,3,2];
$gm_yy += 1900; $gm_mm++;

my %opt;
getopts('de', \%opt) or usage();
my $debugging = $opt{'d'};
my $enhanced_behavior = $opt{'e'};

if (@ARGV) {
  if ( (@ARGV != 1) ||
       (length($ARGV[0]) & 1) || !length($ARGV[0]) || (length($ARGV[0]) > 10)
  ) {
    usage();
  }

  ($hh, $dd, $mm, $yy, $cc) = reverse($ARGV[0] =~ /(..)/g);

  defined($cc) || ($cc = substr($gm_yy,0,2));
  defined($yy) && ($yy = $cc . $yy);
}

if (defined $hh) {
  if (checknum($hh) > 23) {
    warn "$Program: bad hour: '$hh'\n";
    exit EX_FAILURE;
  }
} else {
  $hh = $gm_hh;
}

if (defined $dd) {
  $dd = checknum($dd);
  if ($dd == 0 || $dd > 31) {
    warn "$Program: bad day: '$dd'\n";
    exit EX_FAILURE;
  }
} else {
  $dd = $gm_dd;
}

if (defined $mm) {
  $mm = checknum($mm);
  if ($mm == 0 || $mm > 12) {
    warn "$Program: bad month: '$mm'\n";
    exit EX_FAILURE;
  }
} else {
  $mm = $gm_mm;
}

defined($yy) || ($yy = $gm_yy);

#
## Convert phase time to Julian

sub mdy_to_julian {
  my ($mm, $dd, $yy) = @_;
  my ($j_a, $j_b, $j_c, $j_d);

  my $j_mm = $mm;
  my $j_yy = $yy;

  $j_yy++ if ($j_yy < 0);
  if ($j_mm < 3) {
    $j_mm += 12;
    $j_yy--;
  }
  if ( ($yy < 1582) ||
       ($yy == 1582 and $mm < 10) ||
       ($yy == 1582 and $mm == 10 and $dd < 15)
  ) {
    $j_b = 0;
  }
  else {
    $j_a = floor($j_yy / 100);
    $j_b = 2 - $j_a + floor($j_a / 4);
  }
  if ($j_yy >= 0) {
    $j_c = floor(365.25 * $j_yy) - 694025;
  }
  else {
    $j_c = floor((365.25 * $j_yy) - 0.75) - 694025;
  }
  $j_d = floor(30.6001 * ($j_mm + 1));

  $j_b + $j_c + $j_d + $dd + 2415020;
}

#
## Calculate the eccentric anomaly using Kepler's equation.

sub kepler {
  my ($m, $ecc) = @_;

  my $EPSILON = 1e-6;

  my $e = $m = deg2rad($m);

  my $delta;
  do {
    $delta = $e - $ecc * sin($e) - $m;
    $e -= $delta / (1 - $ecc * cos($e));
  } while (abs($delta) > $EPSILON);

  $e;
}

#
## Calculate the phase of the moon as a fraction (0.0 -> 0.99).  New
## Moon is 0, First Quarter is 0.25, Full Moon is 0.5, Last Quarter is
## 0.75.

sub calc_phase {
  my ($hh, $mm, $dd, $yy) = @_;

  # $utc_offset is in fractions of a day.  For more accurate moon
  # phases, calculate this from the local timezone.

  my $utc_offset = 0;
  my $pdate = mdy_to_julian($mm, $dd, $yy) + ($hh/24) + $utc_offset;

  # sun's position

  my $day = $pdate - EPOCH;
  my $sun_mean_anomaly = fixangle((360 / 365.2422) * $day);
  my $sun_epoch_coords = fixangle($sun_mean_anomaly + ELONGE - ELONGP);

  my $sun_ecc = kepler($sun_epoch_coords, ECCENT);
  $sun_ecc = sqrt((1 + ECCENT) / (1 - ECCENT)) * tan($sun_ecc / 2);
  $sun_ecc = 2 * rad2deg(atan($sun_ecc)); # true anomaly
                                        # sun's geocentric ecliptic longitude
  my $sun_lambda = fixangle($sun_ecc + ELONGP);

  # moon's position

  my $moon_mean_longitude = fixangle(13.1763966 * $day + MMLONG);
  my $moon_mean_anomaly =
    fixangle($moon_mean_longitude - 0.1114041 * $day - MMLONGP);
  my $moon_evection = 1.2739
    * sin(deg2rad(2 * $moon_mean_longitude - $sun_lambda) - $moon_mean_anomaly);
  my $moon_annual_equation = 0.1858 * sin(deg2rad($sun_epoch_coords));
  my $moon_correction_1 = 0.37 * sin(deg2rad($sun_epoch_coords));
  my $moon_corrected_anomaly = $moon_mean_anomaly + $moon_evection -
    $moon_annual_equation - $moon_correction_1;
  my $moon_correction_for_center =
    6.2886 * sin(deg2rad($moon_corrected_anomaly));
  my $moon_correction_2 = 0.214 * sin(deg2rad(2 * $moon_corrected_anomaly));
  my $moon_corrected_longitude =
    $moon_mean_longitude + $moon_evection + $moon_correction_for_center
    - $moon_annual_equation - $moon_correction_2;
  my $moon_variation = 0.6583
    * sin(deg2rad(2 * ($moon_corrected_longitude - $sun_lambda)));
  my $moon_true_longitude = $moon_corrected_longitude + $moon_variation;

  # age of moon, in degrees

  my $moon_age = $moon_true_longitude - $sun_lambda;

  # age of moon as a fraction of a circle

  fixangle($moon_age) / 360;
}

if ($debugging) {
  print "entered    = yy($yy)  mm($mm)  dd($dd)  hh($hh)\n";
  print "julian day = ", mdy_to_julian($mm, $dd, $yy), "\n";
  print "moon phase = ", calc_phase($hh, $mm, $dd, $yy), "\n";
}

#
## Draw the moon.

my $terminal_width = ($ENV{COLS} || $ENV{COLUMNS} || COLS_DEFAULT) - 2;
if (checknum($terminal_width, 'terminal width') == 0) {
  warn "$Program: terminal width 0\n";
  exit EX_FAILURE;
}
my $terminal_height = ($ENV{LINES} || $ENV{ROWS} || LINES_DEFAULT) - 3;
if (checknum($terminal_height, 'terminal height') == 0) {
  warn "$Program: terminal height 0\n";
  exit EX_FAILURE;
}
my @surface_char = (':', '-', '+', '=', '*', 'O', '0', '8', '#', '@');

sub display_moon {
  my ($mm, $dd, $yy, $hh) = @_;

  my $real_phase = calc_phase($hh, $mm, $dd, $yy);
  my $rotated_scaled_phase = ($real_phase * 360 - 180) / 180;
  $real_phase = sprintf('%.2f', $real_phase);
                                        # convert visible phase to percent of full
  my $phase_percent = int($rotated_scaled_phase * 100);
  $phase_percent = ( ($phase_percent < 0)
                     ? (100 + $phase_percent)
                     : (100 - $phase_percent)
                   );
                                        # convert phase percent to name
  my $phase_name;
  if ($phase_percent < 2) {
    $phase_name = 'New';
  }
  elsif ($phase_percent > 99) {
    $phase_name = 'Full';
  }
  elsif ($phase_percent == 50) {
    if ($real_phase < 0.5) {
      $phase_name = 'First Quarter';
    }
    else {
      $phase_name = 'Last Quarter';
    }
  }
  elsif ($phase_percent < 50) {
    if ($real_phase < 0.5) {
      $phase_name = 'Waxing Crescent';
    }
    else {
      $phase_name = 'Waning Crescent';
    }
  }
  else {
    if ($real_phase < 0.5) {
      $phase_name = 'Waxing Gibbous';
    }
    else {
      $phase_name = 'Waning Gibbous';
    }
  }
                                        # happyfun enhanced behavior
  if ($enhanced_behavior) {
    my $increment = 2 / $terminal_height;
    my $count = 0;
    for (my $y=1-($increment/2); $y>-1; $y -= $increment) {
      my $x = sqrt(1 - $y**2);

      my $moon_size  = int($x * $terminal_width * ASPECT_RATIO);
      ($moon_size & 1) && ($moon_size--); # symmetry fudge
      my $space_size = ($terminal_width - $moon_size) / 2;

      $count++;
      my $text = '';
      if ($count==1) {
        $text = "Date : $mm/$dd/$yy:$hh";
      }
      elsif ($count == 2) {
        $text = "Phase: $phase_name";
      }
      elsif ($count == 3) {
        $text = "Full : $phase_percent\%";
      }

      my $space = $text . (' ' x ($space_size - length($text)));

      my $dark_size          = $moon_size * abs($rotated_scaled_phase);
      my $light_size         = $moon_size - $dark_size;
      my $terminus_intensity = $light_size - int($light_size);

      my $dark_side  = $surface_char[0] x $dark_size;
      my $terminus   = $surface_char[$terminus_intensity * @surface_char];
      my $light_side = $surface_char[-1] x $light_size;

      print( $space, ( ($rotated_scaled_phase < 0) ?
                       ($dark_side . $terminus . $light_side) :
                       ($light_side . $terminus . $dark_side)
                     ),
             "\n"
           );
    }
  }
                                        # standard, boring behavior
  else {
    print "The Moon is $phase_name ($phase_percent\% of Full)\n";
  }
}

#
## If debugging, cycle through all the phases for the month.  Compare
## vs. <http://www.googol.com/moon/>.  Otherwise, just display the one
## date/time.

if ($debugging) {
  print SCREEN_CLEAR if ($enhanced_behavior);
  for ($dd = 1; $dd < 32; $dd++) {
    for ($hh = 0; $hh < 24; $hh++) {
      print CURSOR_HOME if ($enhanced_behavior);
      &display_moon($mm, $dd, $yy, $hh);
    }
  }
}
else {
  &display_moon($mm, $dd, $yy, $hh);
}

exit EX_SUCCESS;

__END__

=head1 NAME

  pom - display the phase of the moon

=head1 SYNOPSIS

  pom [-d] [-e] [[[[[[cc]yy]mm]dd]HH]]

=head1 DESCRIPTION

The pom utility displays the current phase of the moon.  This is
useful for selecting software completion target dates, predicting
managerial behavior, or determining thoth's current mood. :)

Pom will display the current moon phase, unless a new time is
specified on the command line.  The parameter's format is similar to
the canonical representation used by date(1).

Pom has two modes: "standard" and "enhanced".  Standard mode presents
a brief description of the moon's phase.  It is the default.  Enhanced
mode shows the phase's date and time, the brief description, and an
ASCII art representation of the moon's phase.

Current options:

  -d    Enable debugging.  This calculates phase information every
        hour for a month.

  -e    Enable enhanced mode.  Turns on the ASCII art feature.

=head1 SEE ALSO

  date(1)

=head1 BUGS

Times must be within the range of the Unix epoch.

The local timezone and coordinates are not considered.

Day validation does not consider the current month, or leap days.

The default character cell aspect ratio is likely to be wrong.

There is some jitter in the current calculations.  It may only be
noticeable while debugging.

The calculations in this version are a few percent different from the
original pom.

=head1 ACKNOWLEDGEMENTS

A lot of the documentation was cribbed from the OpenBSD Reference
Manual.

The moon phase calculations were ported from pcal, which is available
from <ftp://www.decus.org/pub/lib/vs0174/pcal/>.

=head1 AUTHOR and COPYRIGHT

Pom is Copyright 1999 Rocco Caputo <troc@netrus.net>.  All rights
reserved.  Pom is free software; you may redistribute it and/or modify
it under the same terms as Perl itself.