#!/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: Dan Jacobson
# Last Modified On: Thu Aug  7 05:28:08 2003
# Update Count    : 45
 
# 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.

# Invoking this again under perl, if the "#!" line at top didn't do it
# already: [(POSIX) PATHs should have explicit '.'s to make this work
# 100%, i.e., If your PATH ends with "...bla:" you should change it to
# "...bla:."]
eval 'exec perl -S $0 ${1+"$@"}'
    if $running_under_some_shell;

$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.)
";

$start_time = time();

# adjust to 0h UT of nearest day
($gsec, $gmin, $ghr, @rest) = 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);
}
($gsec, $gmin, $ghr, $d, $m, $y, @rest) = gmtime($start_time);
$y += 1900;
$m++;
# Compute jd (Julian date) now
$jd_now = &mdy_to_jd($m,$d,$y);

$stepsize = 1;
$nsteps = 24;
$flip = 1;
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 :
$local_please=1;

# Compute jd
$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;
$[ = 1; # start column numbering at 1
$\ = "\n";

$pi = atan2(0,-1);

#DJ: try localtime too:
($gsec, $gmin, $ghr, $gmd, $gmon, $gyr, @rest) = 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, @rest) = gmtime($start_time);
printf "%02d:%02d:%02d UT on %02d/%02d/%04d\n",
    $ghr, $gmin, $gsec, $gmon+1, $gmd, $gyr + 1900; #Y2K -DJ

($lsec, $lmin, $lhr, @rest) = localtime($start_time); #DJ:

for ($i = 0; $i < $nsteps; $i++) {
# From Chapter 35 and 36 of Meeus
    $d = $jd - 2415020.0 + $i*$stepsize/24;
    $V = 134.63 + 0.00111587 * $d;

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

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

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

    $f1 = $psi - $B;
    $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);
    };

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

    $g = 187.3 + 50.310674 * $f;
    $h = 311.1 + 21.569229 * $f;
  
    $cor_u_1 =  0.472 * &dsin(2*($u_1 - $u_2));
    $cor_u_2 =  1.073 * &dsin(2*($u_2 - $u_3));
    $cor_u_3 =  0.174 * &dsin($g);
    $cor_u_4 =  0.845 * &dsin($h);
  
    $r_1 = 5.9061 - 0.0244 * &dcos(2*($u_1 - $u_2));
    $r_2 = 9.3972 - 0.0889 * &dcos(2*($u_2 - $u_3));
    $r_3 = 14.9894 - 0.0227 * &dcos($g);
    $r_4 = 26.3649 - 0.1944 * &dcos($h);
  
    $x_1 = $r_1 * &dsin($u_1 + $cor_u_1);
    $x_2 = $r_2 * &dsin($u_2 + $cor_u_2);
    $x_3 = $r_3 * &dsin($u_3 + $cor_u_3);
    $x_4 = $r_4 * &dsin($u_4 + $cor_u_4);

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

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

    $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");
$white="        ";

($gsec, $gmin, $ghr, $gmd, $gmon, $gyr, @rest) =
	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, @rest) =
	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 {
    local($month, $day,$year) = @_;
    local($m) = $month;
    local($y) = ($year < 0) ? $year + 1 : $year;

    local($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 {
    local($t) = pop(@_);
    return sin($t*$pi/180);
}

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