#!/bin/sh -- # perl
# \
eval 'exec perl -S $0 ${1+"$@"}'
    if 0;
   
# reason:/data/ast6/1999GU3/analysis/elements.proc

# Parse either email from JPL/Horizons to extract osculating orbital elements 
# as a function of time for a given object.  

# There are two files: one that contains the orbital elements, and the second
# that contains the geocentric range.  Unfortunately there doesn't seem to be
# a way to output everything into one file using Horizons.

# Use Horizons to output the quantities and then email them to yourself.  The
# script below strips the headers (51 lines for a, e, i and 57 lines for delta),
# then it parses for the desired quantities, and finally it outputs the results
# to a separate file for each one (necessary because tkgraph can't handle 
# multiple inputs).

# This software was modified from core code originally written to compute delta-v
# for NEAs using files downloaded from the MPC NEO website. 2000 February 17

# Lance A. M. Benner
# Created:       1999 January 19 (deltavproc)
# elements.proc: 2000 February 17
# Modified:      2000 February 18

print "\n";

# Input files
#$filename[1] = "1999JM8.1979_1988.elements";
#$filename[1] = "1999JM8.1989_1998.elements";
#$filename[1] = "1999JM8.1999_2008.elements";
$filename[1] = "1999JM8.2009_2018.elements";

$filename[2] = "1999JM8.1979_2018.range";

$filemin = 1;
$filemax = 1;

print "Reading the input files:\n";
for ($file = $filemin; $file <= $filemax; $file++)
{
	print "  $filename[$file]\n";	# Verify that input files are being read
	@lines = ();

# Parse the input files and output the results  
	open(FILE, $filename[$file]) || die "Can't open $filename[$file]: $!\n";
	@lines = <FILE>;
	close(FILE);
#	print "@lines";                 # Verify that all input lines were read
	$j = 0;
	$n = 0;
	foreach (@lines)                # Parse each input line
	{
# Parse lines past the header counting from zero
	  if($file == 1)                 # Select for a, e, i
	    {	      
		if($j > 50 && /(.{17}).{33}(.{23}).{1}(.{23}).{1}(.{23}).+/)	# Ignore the header
		{
			$date[$n] = $1;			# Julian date
			$e[$n] = $2;			# Eccentricity
			$qr[$n] = $3;			# Perihelion (AU)
			$in[$n] = $4;			# Inclination (deg)
			$a[$n] = $qr[$n]/(1.0 - $e[$n]);# Semimajor axis (AU)
			$n++;	# Number of lines of data in all files so far
		}
	    }

	  if($file == 2)                 # Select for the range (delta)
	    {	      
		if($j > 56 && /(.{16})\s+(.{12}).+/)	# Ignore the header
		{
			$date[$n] = $1;			# Julian date
			$delta[$n] = $2;		# Delta (AU)
#                        if($j < 100)
#                        {
#			     print "$date[$n]  $delta[$n]\n";
#        		}
			$n++;	# Number of lines of data in all files so far
		}
	    }
	    $j++; 			# Number of lines per file so far
	} 				# End of @lines loop
	print "\n";
} 					# End of the files loop

# Output the results
print "  \n";
print "...writing the results to a file...\n";

open(TABLE1,">semimajoraxis.dat") || die "Can't create";
for ($k = 1; $k < $n; $k++)
{
	write TABLE1;
}
close(TABLE1);

open(TABLE2,">eccentricity.dat") || die "Can't create";
for ($k = 1; $k < $n; $k++)
{
	write TABLE2;
}
close(TABLE2);

open(TABLE3,">inclination.dat") || die "Can't create";
for ($k = 1; $k < $n; $k++)
{
	write TABLE3;
}
close(TABLE3);

open(TABLE4,">delta.dat") || die "Can't create";
for ($k = 1; $k < $n; $k++)
{
	write TABLE4;
}
close(TABLE4);

print "...done\n";
print "\n";

########################################################################################
#
#                               Formats
#
########################################################################################

format TABLE1 = 
@#######.##     @##.##### 
$date[$k], $a[$k]
.

format TABLE2 = 
@#######.##     @#.######
$date[$k], $e[$k]
.

format TABLE3 = 
@#######.##     @##.######
$date[$k], $in[$k]
.

format TABLE4 = 
@#######.##     @##.######
$date[$k], $delta[$k]
.

# --------------------------------------------------------------------