[sci.astro] Jupiter's major moons simple plotter in perl

ccount@ATHENA.MIT.EDU, Dan_Jacobson@ATT.COM (05/02/90)

Here's a  perl program to   calculate the positions  of  the Galilean
satellites of Jupiter, and display them in 80 columns of text.  It was
inspired by Dan's awk  version posted briefly  in March, but uses more
accurate formulas     from   Meeus's  "Astronomical      Formulae  for
Calculators".

This is a beta+ test, before submission to comp.sources.misc.
Improvements and suggestions are sought.

Here's an example:
jupmoons nsteps=13 stepsize=2 d=3 reverse
Julian day: 2448014.500000
 0:00:00 UT on 5/3/90
                      g               J00     i  e            c
                        g             J02      i  e          c
                         g            J04      i   e        c
                          g           J06      i    e      c
                            g         J08     i      e    c
                             g        J10   i        e   c
                               g      J12 i          e  c
                                g     Ji4           e  c
                                  g  iJ16          e c
                                  ig  J18         e c
                                i    gJ20        e c
                               i      Jg2       e c
                              i       J0g     e  c
W                            0:00:00 UT on  5/ 4/90                           E


Enjoy,
Craig Counterman and Dan Jacobson

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  jupmoons
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'jupmoons' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'jupmoons'\"
else
echo shar: Extracting \"'jupmoons'\" \(7565 characters\)
sed "s/^X//" >'jupmoons' <<'END_OF_FILE'
X#!/usr/bin/perl
X# jupmoons, version 5/1/1990
X
X# This is jupmoons: a poorperson's Jupiter's main moons simple ASCII
X# grapher for UNIX.  Perl version by Craig Counterman
X# (ccount@athena.mit.edu), inspired by version by Dan_Jacobson@ATT.COM
X# 3/1990 (who helped a bit on this version too).  Formulae from
X# "Astronomical Formulae for Calculators", by Jean Meeus.  Richmond,
X# Va., U.S.A. : Willmann-Bell, c1982.
X
X# Copyright (c) 1990 by Craig Counterman and Dan Jacobson. All rights
X# reserved. See the end of this file for details.
X
X# Make simple ASCII graph of Jupiter and the relative positions of its
X# 4 major moons.  Each moon is represented by its initial letter,
X# Jupiter is located at the three digit number (hour) in the center.
X# The scale is 1.5 char per Jupiter radius, center of Jupiter is
X# column 40.  A usage description is below.
X
X# Possible additional feature ideas: output (and input?) in local
X# time.
X
X# Invoking this again under perl, if the "#!" line at top didn't do it
X# already: [(POSIX) PATHs should have explicit '.'s to make this work
X# 100%, i.e., If your PATH ends with "...bla:" you should change it to
X# "...bla:."]
Xeval 'exec perl -S $0 ${1+"$@"}'
X    if $running_under_some_shell;
X
X$USAGE = "$0: usage:
XArguments M=m, D=d, Y=yr, default to current GMT date + 1 day.
XSTEPSIZE=hrs in hours, default 1, NSTEPS=n, default 24. Specify
XREVERSE to make left west.  (Arguments can be in upper or lower case.)
X";
X
X$start_time = time();
X
X# adjust to 0h UT of nearest day
X($gsec, $gmin, $ghr, @rest) = gmtime($start_time);
Xif (($ghr*60*60 + $gmin*60 + $gsec) > 12*60*60) {
X    $start_time +=  24*60*60 - ($ghr*60*60 + $gmin*60 + $gsec);
X} else {
X    $start_time -=  ($ghr*60*60 + $gmin*60 + $gsec);
X}
X($gsec, $gmin, $ghr, $d, $m, $y, @rest) = gmtime($start_time);
X$y += 1900;
X$m++;
X# Compute jd (Julian date) now
X$jd_now = &mdy_to_jd($m,$d,$y);
X
X$stepsize = 1;
X$nsteps = 24;
X$flip = 1;
Xfor (@ARGV) {
X    if (/^stepsize=[\d.]+$/i) {
X	@str = /^stepsize= *([\d.]*)/i;
X	$stepsize = @str[0];
X    } elsif (/^nstep[s]?=\d+$/i) {
X# we're accepting nstep or nsteps, though nsteps is documented, 
X# because Craig found he typed nstep sometimes by mistake
X	@str = /^nstep[s]?= *(\d*)/i;
X	$nsteps = @str[0];
X    } elsif (/^m=\d+$/i) {
X	@str = /^m= *(\d*)/i;
X	$m = @str[0];
X    } elsif (/^d=\d+$/i) {
X	@str = /^d= *(\d*)/i;
X	$d = @str[0];
X    } elsif (/^y=\d+$/i) {
X	@str = /^y= *(\d*)/i;
X	$y = @str[0];
X    } elsif ((/^reverse$/i) || (/^we$/i)) {
X	$flip = -1;
X    } else {
X	die $USAGE;
X    }
X}
X
X# Compute jd
X$jd = &mdy_to_jd($m,$d,$y);
X# compute revised $start_time
X$start_time += ($jd - $jd_now)*86400;
X
Xprintf "Julian day: %f\n", $jd;
X$[ = 1; # start column numbering at 1
X$\ = "\n";
X
X$pi = atan2(0,-1);
X
X
X($gsec, $gmin, $ghr, $gmd, $gmon, $gyr, @rest) = gmtime($start_time);
Xprintf "%2d:%02d:%02d UT on %d/%d/%d\n",
X    $ghr, $gmin, $gsec, $gmon+1, $gmd, $gyr;
X
Xfor ($i = 0; $i < $nsteps; $i++) {
X# From Chapter 35 and 36 of Meeus
X    $d = $jd - 2415020.0 + $i*$stepsize/24;
X    $V = 134.63 + 0.00111587 * $d;
X
X    $M = (358.47583 + 0.98560003*$d);
X    $N = (225.32833 + 0.0830853*$d) + 0.33 * &dsin($V);
X
X    $J = 221.647 + 0.9025179*$d - 0.33 * &dsin($V);;
X
X    $A = 1.916*&dsin($M)+0.02*&dsin(2*$M);
X    $B =5.552*&dsin($N)+0.167*&dsin(2*$N);
X    $K = ($J+$A-$B);
X    $R = 1.00014 - 0.01672 * &dcos($M) - 0.00014 * &dcos(2*$M);
X    $r = 5.20867 - 0.25192 * &dcos($N) - 0.00610 * &dcos(2*$N);
X    $Delta = sqrt($R*$R + $r*$r - 2*$R*$r*&dcos($K));
X    $z = $R/$Delta*&dsin($K);
X    $psi = atan2($z, sqrt(1-$z*$z))*180/$pi;
X    $f = ($d - $Delta*0.00577167643528);	# speed of light correction
X
X    $f1 = $psi - $B;
X    $outline = " " x 80;
X
X# Jupiter
X    substr($outline, 39, 3) =  "J" . sprintf("%02.0f",
X                                       ($ghr + $gmin/60 + $i*$stepsize) % 24);
X
X    $u_1 = 84.5506 + 203.4058630 * $f + $f1;
X    $u_2 = 41.5015 + 101.2916323 * $f + $f1;
X    $u_3 = 109.9770 + 50.2345169 * $f + $f1;
X    $u_4 = 176.3586 + 21.4879802 * $f + $f1;
X
X    $g = 187.3 + 50.310674 * $f;
X    $h = 311.1 + 21.569229 * $f;
X  
X    $cor_u_1 =  0.472 * &dsin(2*($u_1 - $u_2));
X    $cor_u_2 =  1.073 * &dsin(2*($u_2 - $u_3));
X    $cor_u_3 =  0.174 * &dsin($g);
X    $cor_u_4 =  0.845 * &dsin($h);
X  
X    $r_1 = 5.9061 - 0.0244 * &dcos(2*($u_1 - $u_2));
X    $r_2 = 9.3972 - 0.0889 * &dcos(2*($u_2 - $u_3));
X    $r_3 = 14.9894 - 0.0227 * &dcos($g);
X    $r_4 = 26.3649 - 0.1944 * &dcos($h);
X  
X    $x_1 = $r_1 * &dsin($u_1 + $cor_u_1);
X    $x_2 = $r_2 * &dsin($u_2 + $cor_u_2);
X    $x_3 = $r_3 * &dsin($u_3 + $cor_u_3);
X    $x_4 = $r_4 * &dsin($u_4 + $cor_u_4);
X
X    $z_1 = $r_1 * &dcos($u_1 + $cor_u_1);
X    $z_2 = $r_2 * &dcos($u_2 + $cor_u_2);
X    $z_3 = $r_3 * &dcos($u_3 + $cor_u_3);
X    $z_4 = $r_4 * &dcos($u_4 + $cor_u_4);
X
X    $lambda = 238.05 + 0.083091*$d + 0.33 * &dsin($V) + $B;
X    $D_s = 3.07*&dsin($lambda+44.5);
X    $D_e = $D_s - 2.15 * &dsin($psi) * &dcos($lambda+24) 
X	- 1.31 * ($r - $Delta)*&dsin($lambda - 99.4)/$Delta;
X    $y_1 = - $r_1 * &dcos($u_1 + $cor_u_1) * &dsin($D_e);
X    $y_2 = - $r_2 * &dcos($u_2 + $cor_u_2) * &dsin($D_e);
X    $y_3 = - $r_3 * &dcos($u_3 + $cor_u_3) * &dsin($D_e);
X    $y_4 = - $r_4 * &dcos($u_4 + $cor_u_4) * &dsin($D_e);
X
X    $t_4 = 40 +$flip * 1.5 * $x_4;
X    if ($t_4 < 1) {$t_4 = 1;}
X    substr($outline, (40 +$flip * 1.5 * $x_1), 1) = "i";
X    substr($outline, (40 +$flip * 1.5 * $x_2), 1) = "e";
X    substr($outline, (40 +$flip * 1.5 * $x_3), 1) = "g";
X    substr($outline, $t_4, 1) = "c";
X
X    $outline =~ s/ *$//;
X    print $outline;
X}
X
X
X($gsec, $gmin, $ghr, $gmd, $gmon, $gyr, @rest) =
X	gmtime($start_time+($i-1)*$stepsize*60*60);
Xif ($flip == 1) {
X    printf "E%s%2d:%02d:%02d UT on %2d/%2d/%2d%sW\n", " " x 27,
X    $ghr, $gmin, $gsec, $gmon+1, $gmd, $gyr, " " x 27;
X} else {
X    printf "W%s%2d:%02d:%02d UT on %2d/%2d/%2d%sE\n", " " x 27,
X    $ghr, $gmin, $gsec, $gmon+1, $gmd, $gyr, " " x 27;
X}
X
Xsub mdy_to_jd {
X    local($month, $day,$year) = @_;
X    local($m) = $month;
X    local($y) = ($year < 0) ? $year + 1 : $year;
X
X    local($a, $b, $c, $d);
X    $, = " ";
X    $\ = "\n";
X
X    if ($month < 3) {
X	$m += 12;
X	$y -= 1;
X    }
X
X    if (($year < 1582) ||
X	($year == 1582 && (($month < 10)
X			   || (($month == 10) && ($day < 15))))) {
X	$b = 0;
X    } else {
X	$a = int($y/100);
X	$b = 2 - $a + int($a/4);
X    }
X
X    if ($y < 0) {
X	$c = int((365.25*$y) - 0.75) + 1720995;
X    } else {
X	$c = int(365.25*$y) + 1720995;
X    }
X
X    $d = int(30.6001*($m+1));
X
X    return ($b + $c + $d + $day - 0.5);
X}
X
Xsub dsin {
X    local($t) = pop(@_);
X    return sin($t*$pi/180);
X}
X
Xsub dcos {
X    local($t) = pop(@_);
X    return cos($t*$pi/180);
X}
X
X# Copyright (c) 1990 by Craig Counterman and Dan Jacobson. All rights
X# reserved
X
X# This software may be redistributed freely, not sold.  The copyright
X# notices and disclaimer of warranty must remain unchanged.
X
X# No representation is made about the suitability of this software for
X# any purpose.  It is provided "as is" without express or implied
X# warranty, to the extent permitted by applicable law.
X
X# DISCLAIMER OF WARRANTY
X# ----------------------
X# The authors disclaim all warranties with regard to this software to
X# the extent permitted by applicable law, including all implied
X# warranties of merchantability and fitness. In no event shall the
X# authors be liable for any special, indirect or consequential damages
X# or any damages whatsoever resulting from loss of use, data or
X# profits, whether in an action of contract, negligence or other
X# tortious action, arising out of or in connection with the use or
X# performance of this software.
END_OF_FILE
if test 7565 -ne `wc -c <'jupmoons'`; then
    echo shar: \"'jupmoons'\" unpacked with wrong size!
fi
chmod +x 'jupmoons'
# end of 'jupmoons'
fi
echo shar: End of shell archive.
exit 0