#!/usr/bin/perl
# jupmoons -- simple ASCII Jupiter moon grapher
# Author          : Craig Counterman (ccount@mit.edu)
# Author          : Dan Jacobson http://jidanni.org/
# Copyright       : http://www.fsf.org/copyleft/gpl.html
# Created On      : 1990
# Last Modified By: jidanni
# Last Modified On: Fri Dec 13 07:51:12 2013
# Update Count    : 71

# 2001 Nov: Dan Jacobson making minor adjustments to the
# ftp://ftp.uu.net/usenet/comp.sources.misc/volume13/jupmoons
# version of 1990.

# This is jupmoons: a poorperson's Jupiter's main moons simple ASCII
# grapher.  Perl version May 1 1990 by Craig Counterman
# (ccount@mit.edu), inspired by version by Dan Jacobson 3/1990 (who
# helped a bit on this version too).  Formulae from "Astronomical
# Formulae for Calculators", by Jean Meeus.  Richmond, Va., U.S.A. :
# Willmann-Bell, c1982.

# Copyright (c) 1990 by Craig Counterman and Dan Jacobson.
# Copyright (c) 2001 FSF, see http://www.fsf.org/copyleft/gpl.html

# Make simple ASCII graph of Jupiter and the relative positions of its
# 4 major moons.  Each moon is represented by its initial letter,
# Jupiter is located at the three digit number (hour) in the center.
# The scale is 1.5 char per Jupiter radius, center of Jupiter is
# column 40.  A usage description is below.

# Possible additional feature ideas: output (and input?) in local
# time.

# Dan, 2001, says: yes indeed, we simple ascii viewer folks deserve
# local time as the default!  I tried to add some stuff, but not much,
# as I just can't hack it. --- OK, hacked it in... but didn't make it
# controllable from the command line because of weak perl skills.
# Users within 90 degrees of longitude of Greenwich still find their
# night's display cut in half with the default settings...
# Anyway, for now I can only call it a "quasi"-night-at-a-glance lister.

# really ought to do: make it list from 5PM-7AM local time by default...
# or even exclude times when jupiter is not up... but thats getting too
# complicated ... ok end of 2001 comments.

use strict;
use warnings FATAL => 'all';

my $USAGE = "$0: usage:
Arguments M=m, D=d, Y=yr, default to current GMT date + 1 day.
STEPSIZE=hrs in hours, default 1, NSTEPS=n, default 24. Specify
REVERSE to make left west.  (Arguments can be in upper or lower case.)
";

my $start_time = time();

# adjust to 0h UT of nearest day
my ( $gsec, $gmin, $ghr, undef ) = gmtime($start_time);
if ( ( $ghr * 60 * 60 + $gmin * 60 + $gsec ) > 12 * 60 * 60 ) {
    $start_time += 24 * 60 * 60 - ( $ghr * 60 * 60 + $gmin * 60 + $gsec );
}
else {
    $start_time -= ( $ghr * 60 * 60 + $gmin * 60 + $gsec );
}
my ( $d, $m, $y );
( $gsec, $gmin, $ghr, $d, $m, $y, undef ) = gmtime($start_time);
$y += 1900;
$m++;

# Compute jd (Julian date) now
my $jd_now = &mdy_to_jd( $m, $d, $y );

my $stepsize = 1;
my $nsteps   = 24;
my $flip     = 1;
my @str;
for (@ARGV) {
    if (/^stepsize=[\d.]+$/i) {
        @str      = /^stepsize= *([\d.]*)/i;
        $stepsize = $str[0];
    }
    elsif (/^nstep[s]?=\d+$/i) {

        # we're accepting nstep or nsteps, though nsteps is documented,
        # because Craig found he typed nstep sometimes by mistake
        @str    = /^nstep[s]?= *(\d*)/i;
        $nsteps = $str[0];
    }
    elsif (/^m=\d+$/i) {
        @str = /^m= *(\d*)/i;
        $m   = $str[0];
    }
    elsif (/^d=\d+$/i) {
        @str = /^d= *(\d*)/i;
        $d   = $str[0];
    }
    elsif (/^y=\d+$/i) {
        @str = /^y= *(\d*)/i;
        $y   = $str[0];
    }
    elsif ( (/^reverse$/i) || (/^we$/i) ) {
        $flip = -1;
    }
    else {
        die $USAGE;
    }
}

#DJ: perhaps future command line [L=0 to turn off] switch, hardwired for now :
my $local_please = 1;

# Compute jd
my $jd = &mdy_to_jd( $m, $d, $y );

# compute revised $start_time
$start_time += ( $jd - $jd_now ) * 86400;

printf "Jupiter's moons.  Julian day: %f\n", $jd;

my $pi = atan2( 0, -1 );

#DJ: try localtime too:
my ( $gmd, $gmon, $gyr );
( $gsec, $gmin, $ghr, $gmd, $gmon, $gyr, undef ) = localtime($start_time);
printf "         %02d:%02d:%02d local on %02d/%02d/%02d == ",
  $ghr, $gmin, $gsec, $gmon + 1, $gmd, $gyr + 1900;

( $gsec, $gmin, $ghr, $gmd, $gmon, $gyr, undef ) = gmtime($start_time);
printf "%02d:%02d:%02d UT on %02d/%02d/%04d\n",
  $ghr, $gmin, $gsec, $gmon + 1, $gmd, $gyr + 1900;    #Y2K -DJ

my ( $lsec, $lmin, $lhr, undef ) = localtime($start_time);    #DJ:
my $i;
for ( $i = 0 ; $i < $nsteps ; $i++ ) {

    # From Chapter 35 and 36 of Meeus
    my $d = $jd - 2415020.0 + $i * $stepsize / 24;
    my $V = 134.63 + 0.00111587 * $d;

    my $M = ( 358.47583 + 0.98560003 * $d );
    my $N = ( 225.32833 + 0.0830853 * $d ) + 0.33 * &dsin($V);

    my $J = 221.647 + 0.9025179 * $d - 0.33 * &dsin($V);

    my $A     = 1.916 * &dsin($M) + 0.02 * &dsin( 2 * $M );
    my $B     = 5.552 * &dsin($N) + 0.167 * &dsin( 2 * $N );
    my $K     = ( $J + $A - $B );
    my $R     = 1.00014 - 0.01672 * &dcos($M) - 0.00014 * &dcos( 2 * $M );
    my $r     = 5.20867 - 0.25192 * &dcos($N) - 0.00610 * &dcos( 2 * $N );
    my $Delta = sqrt( $R * $R + $r * $r - 2 * $R * $r * &dcos($K) );
    my $z     = $R / $Delta * &dsin($K);
    my $psi   = atan2( $z, sqrt( 1 - $z * $z ) ) * 180 / $pi;
    my $f = ( $d - $Delta * 0.00577167643528 );    # speed of light correction

    my $f1      = $psi - $B;
    my $outline = " " x 80;

    # Jupiter
    if ( !$local_please ) {                        #DJ
        substr( $outline, 39, 3 ) = "J"
          . sprintf( "%02.0f", ( $ghr + $gmin / 60 + $i * $stepsize ) % 24 );
    }
    else {
        substr( $outline, 39, 3 ) =
          sprintf( "%02.0f", ( $lhr + $lmin / 60 + $i * $stepsize ) % 24 );
    }

    my $u_1 = 84.5506 + 203.4058630 * $f + $f1;
    my $u_2 = 41.5015 + 101.2916323 * $f + $f1;
    my $u_3 = 109.9770 + 50.2345169 * $f + $f1;
    my $u_4 = 176.3586 + 21.4879802 * $f + $f1;

    my $g = 187.3 + 50.310674 * $f;
    my $h = 311.1 + 21.569229 * $f;

    my $cor_u_1 = 0.472 * &dsin( 2 * ( $u_1 - $u_2 ) );
    my $cor_u_2 = 1.073 * &dsin( 2 * ( $u_2 - $u_3 ) );
    my $cor_u_3 = 0.174 * &dsin($g);
    my $cor_u_4 = 0.845 * &dsin($h);

    my $r_1 = 5.9061 - 0.0244 * &dcos( 2 * ( $u_1 - $u_2 ) );
    my $r_2 = 9.3972 - 0.0889 * &dcos( 2 * ( $u_2 - $u_3 ) );
    my $r_3 = 14.9894 - 0.0227 * &dcos($g);
    my $r_4 = 26.3649 - 0.1944 * &dcos($h);

    my $x_1 = $r_1 * &dsin( $u_1 + $cor_u_1 );
    my $x_2 = $r_2 * &dsin( $u_2 + $cor_u_2 );
    my $x_3 = $r_3 * &dsin( $u_3 + $cor_u_3 );
    my $x_4 = $r_4 * &dsin( $u_4 + $cor_u_4 );

    my $z_1 = $r_1 * &dcos( $u_1 + $cor_u_1 );
    my $z_2 = $r_2 * &dcos( $u_2 + $cor_u_2 );
    my $z_3 = $r_3 * &dcos( $u_3 + $cor_u_3 );
    my $z_4 = $r_4 * &dcos( $u_4 + $cor_u_4 );

    my $lambda = 238.05 + 0.083091 * $d + 0.33 * &dsin($V) + $B;
    my $D_s    = 3.07 * &dsin( $lambda + 44.5 );
    my $D_e =
      $D_s -
      2.15 * &dsin($psi) * &dcos( $lambda + 24 ) -
      1.31 * ( $r - $Delta ) * &dsin( $lambda - 99.4 ) / $Delta;
    my $y_1 = -$r_1 * &dcos( $u_1 + $cor_u_1 ) * &dsin($D_e);
    my $y_2 = -$r_2 * &dcos( $u_2 + $cor_u_2 ) * &dsin($D_e);
    my $y_3 = -$r_3 * &dcos( $u_3 + $cor_u_3 ) * &dsin($D_e);
    my $y_4 = -$r_4 * &dcos( $u_4 + $cor_u_4 ) * &dsin($D_e);

    my $t_4 = 40 + $flip * 1.5 * $x_4;
    if ( $t_4 < 1 ) { $t_4 = 1; }
    substr( $outline, ( 40 + $flip * 1.5 * $x_1 ), 1 ) = "i";
    substr( $outline, ( 40 + $flip * 1.5 * $x_2 ), 1 ) = "e";
    substr( $outline, ( 40 + $flip * 1.5 * $x_3 ), 1 ) = "g";
    substr( $outline, $t_4, 1 ) = "c";

    $outline =~ s/ *$//;
    print $outline;
}

#DJ
printf( ( $flip == 1 ) ? "E" : "W" );
my $white = "        ";

( $gsec, $gmin, $ghr, $gmd, $gmon, $gyr, undef ) =
  localtime( $start_time + ( $i - 1 ) * $stepsize * 60 * 60 );
printf "%s%02d:%02d:%02d local on %02d/%02d/%04d == ", $white,
  $ghr, $gmin, $gsec, $gmon + 1, $gmd, $gyr + 1900;

( $gsec, $gmin, $ghr, $gmd, $gmon, $gyr, undef ) =
  gmtime( $start_time + ( $i - 1 ) * $stepsize * 60 * 60 );
printf "%02d:%02d:%02d UT on %02d/%02d/%04d%s",
  $ghr, $gmin, $gsec, $gmon + 1, $gmd, $gyr + 1900, $white;

printf( ( $flip != 1 ) ? "E" : "W" );
print "";

#end of DJ

sub mdy_to_jd {

    my ( $month, $day, $year ) = @_;
    my ($m) = $month;
    my ($y) = ( $year < 0 ) ? $year + 1 : $year;

    my ( $a, $b, $c, $d );
    $, = " ";
    $\ = "\n";

    if ( $month < 3 ) {
        $m += 12;
        $y -= 1;
    }

    if (
        ( $year < 1582 )
        || (
            $year == 1582
            && (   ( $month < 10 )
                || ( ( $month == 10 ) && ( $day < 15 ) ) )
        )
      )
    {
        $b = 0;
    }
    else {
        $a = int( $y / 100 );
        $b = 2 - $a + int( $a / 4 );
    }

    if ( $y < 0 ) {
        $c = int( ( 365.25 * $y ) - 0.75 ) + 1720995;
    }
    else {
        $c = int( 365.25 * $y ) + 1720995;
    }

    $d = int( 30.6001 * ( $m + 1 ) );

    return ( $b + $c + $d + $day - 0.5 );
}

sub dsin {
    my ($t) = pop(@_);
    return sin( $t * $pi / 180 );
}

sub dcos {
    my ($t) = pop(@_);
    return cos( $t * $pi / 180 );
}