Perl script to convert JPL Horizons Ephemeris data to Celestia SSC Elliptical Orbit Parameters

With Comet ISON coming up, I've been playing with Celestia a bit to see if I can get it to show me with reasonable accuracy where the spacecraft SDO and STEREO (-A and -B) are in relation to the comet.

Long story short -- I had to make my own Celestia ".ssc" files to describe the spacecraft. Using the JPL HORIZONS web interfaceyou can get recent orbital elements, but you have to convert these to something Celestia will understand (an "EllipticalOrbit" stanza in the ".ssc" file).

The perl script below does the conversion; just run it and paste in a single block of Ephemeris data like this:

2456621.500000000 = A.D. 2013-Nov-25 00:00:00.0000 (CT)
 EC= 4.143129880700237E-02 QR= 9.996672652291025E-01 IN= 2.936558964716465E-01
 OM= 3.363978705358475E+02 W = 1.469845587993808E+02 Tp=  2456458.185749091208
 N = 9.254559117051611E-01 MA= 1.511401389693867E+02 TA= 1.533315740323342E+02
 A = 1.042874927988944E+00 AD= 1.086082590748784E+00 PR= 3.889974610856358E+02

And you'll get output line this:

EllipticalOrbit
{
        Epoch 2456621.500000000
        Period 1.06499689086802
        SemiMajorAxis 1.04287492798894
        Eccentricity 0.041431298807
        Inclination 0.293655896472
        AscendingNode 336.397870535848
        ArgOfPericenter 146.984558799381
        MeanAnomaly 151.140138969387
}

The perl script is below, or you can download it here.

#!/usr/bin/perl -w

# input an ephem set like this:
#2456621.500000000 = A.D. 2013-Nov-25 00:00:00.0000 (CT)
# EC= 4.143129880700237E-02 QR= 9.996672652291025E-01 IN= 2.936558964716465E-01
# OM= 3.363978705358475E+02 W = 1.469845587993808E+02 Tp=  2456458.185749091208
# N = 9.254559117051611E-01 MA= 1.511401389693867E+02 TA= 1.533315740323342E+02
# A = 1.042874927988944E+00 AD= 1.086082590748784E+00 PR= 3.889974610856358E+02

# output like this:
#        EllipticalOrbit
#        {
#                Epoch 2456621.500000000 # = A.D. 2013-Nov-25 00:00:00.0000 (CT)
#                # reference http://www.lepp.cornell.edu/~seb/celestia/transforming_ephemeris.html
#                Period 1.06499689086801936 # P = (q/(1-e))^1.5 for a closed, elliptical orbit (q = QR, e = ER)
#                SemiMajorAxis 1.04287492798894325 # a = (q+Q)/2, where Q is Apoapsis distance (AD in horizons ELEMENTS)
#                Eccentricity 0.04143129880700237 # (EC)
#                Inclination 0.2936558964716465 # (IN)
#                AscendingNode 336.3978705358475 # (OM)
#                ArgOfPericenter 146.9845587993808 # (W)
#                MeanAnomaly 151.1401389693867 # (MA)
#        }

my ($Epoch, $EC, $QR, $IN, $OM, $W, $Tp, $N, $MA, $TA, $A, $AD, $PR);
my ($Period, $SemiMajorAxis);


print "Input a single block of ephem data\n";

while () {
        if ( /^((\w|\.)+)/ ) {
                $Epoch = $1;
        }

        if ( /EC=/ ) {
                ($EC, $QR, $IN) = (split())[1,3,5];
        }

        if ( /OM= / ) {
                #$_ = readline; chomp;
                ($OM, $W, $Tp) = (split())[1,4,6];
        }

        if ( /N = / ) {
                #$_ = readline; chomp;
                ($N, $MA, $TA) = (split())[2,4,6];
        }

        if ( /AD= / ) {
                #$_ = readline; chomp;
                ($A, $AD, $PR) = (split())[2,4,6];
                last;
        }
}

$Period = ( $QR / (1 - $EC) )**1.5 ;
$SemiMajorAxis = ( $QR + $AD ) / 2 ;

print "EllipticalOrbit\n";
print "{\n";
print "\tEpoch $Epoch\n";
print "\tPeriod $Period\n";
print "\tSemiMajorAxis $SemiMajorAxis\n";
print "\tEccentricity " . sprintf("%12.12f\n",$EC);
print "\tInclination " . sprintf("%12.12f\n",$IN);
print "\tAscendingNode " . sprintf("%12.12f\n",$OM);
print "\tArgOfPericenter " . sprintf("%12.12f\n",$W);
print "\tMeanAnomaly " . sprintf("%12.12f\n",$MA);
print "}\n";