Comets and Minor Planets

[Bernmeister]

This discussion serves two objectives: firstly, to let others know of the changes at the Minor Planet Center (MPC) in terms of orbital element data for comets and minor planets; secondly a cry for help as it were for my own specific case...

Background

I have a Python program which pulls down data files from the MPC for comets and minor planets. My program uses PyEphem and so was grabbing the XEphem specific data files, located here. For each of the data files, I iterate over each body and filter out those with an apparent magnitude greater than 15. This takes ten to twenty seconds.

A few months ago I noticed some discrepancies between the XEphem data and the native MPC data and brought that to the attention of the MPC. Related or not, the XEphem data page now only contains a link for comets; the data files for bright/critical/distant/unusual minor planets are no longer available. In reality, the original URLs are still active, but the hyperlinks are no longer present on the page. The equivalent MPC data page too only contains a link for comet data (all other minor planet data hyperlinks have been stripped but still exist as URLs).

In parallel, I have converted my program to optionally use Skyfield. Whilst all the functionality is there, when I do the corresponding iteration over each minor planet or comet to filter by magnitude, the process can take several minutes under Skyfield. So I am somewhat "stuck" with PyEphem for now.

The Problem

When I brought to the attention of the MPC the discrepancies between XEphem and MPC format data file, I was told that the XEphem files have not likely been updated for a few years. Given that and that only the comet data is now available for XEphem, what to do next...

Solution 1

I could continue to use PyEphem and drop the minor planet functionality.

Solution 2

I could switch to Skyfield and only support comets (there's not too many comets compared to minor planets so the filtering out by magnitude should not be too time consuming). There is an open issue to increase the speed in Skyfield for processing comets and minor planets, but this is a non-trivial task.

Solution 3

Continue with PyEphem using the comet data as still available. For minor planets, although the data files for bright/critical/distant/unusual are no longer available, there are other files that are similar enough but in MPC format. I found on the XEphem page a script to convert MPC data to XEphem format for minor planets. Assuming this program works as I hope it does, I should be able to feed in various MPC format minor planet data files and get XEphem out.

@brandon-rhodes: I know you are not the author of mpcorb2edb.pl, but do you know if this works as intended? If so, I might even be tempted to write a Python version.

Suggestions welcome!

I will post a link to this discussion on the PyEphem discussion page as it is related.

[Bernmeister]

Update

Given that the data files for bright/critical/distant/unusual minor planets are no longer (officially) available and likely have not been updated for some time, I thought about swapping over to use the files on the main data page, namely NEA.txt, PHA.txt, Distant.txt and Unusual.txt. Unfortunately many bodies, such as Ceres, Juno and Vesta are not present in these new files (but are available in the older originals).

The next option was to use the full list of minor planets, namely MPCORB.DAT. I wrote a Python script to read in each line (one body per line), instantiate a body and determine the current apparent magnitude (using PyEphem). Any body which exceeds apparent magnitude of 15 is discarded. As of May 2022, there are around 1.2 million bodies in MPCORB.DAT and after filtering, I end up with around 1200 and a file size of around 250kb, taking around five minutes to complete.

For anyone who needs data for minor planets and wants to screen out bodies by apparent magnitude ("visible" bodies), you are welcome to use the script below. The script reads in MPC formatted data for minor planets, converts to XEphem/PyEphem format, then instantiates each body and filters by apparent magnitude.

Filter MPCORB.DAT by apparent magnitude:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Read in MPCORB.DAT and for each body, compute the apparent magnitude over a period using PyEphem.
# If a body's apparent magnitude is less than the maximum specified at any point, keep the body.
#
# The MCPORB.DAT data which is in MPC format, must be converted to XEphem format.
# Refer to https://github.com/XEphem/XEphem/blob/main/GUI/xephem/tools/mpccomet2edb.pl
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import datetime, ephem


# https://www.minorplanetcenter.net/iau/info/PackedDates.html
centuryMap = { 'I': 1800, 'J': 1900, 'K': 2000 }


# https://www.minorplanetcenter.net/iau/info/PackedDates.html
def getUnpackedDate( i ): return str( i ) if i.isdigit() else str( ord( i ) - ord( 'A' ) + 10 )


def mpcToXEphem( line ):
    xEphemLine = ""
    name = line[ 167 - 1 : 194 ].replace( '(', '' ).replace( ')', '' ).strip()
    absoluteMagnitude = line[ 9 - 1 : 13 ].strip()
    inclinationToEcliptic = line[ 60 - 1 : 68 ].strip()
    longitudeAscendingNode = line[ 49 - 1 : 57 ].strip()
    argumentPerihelion = line[ 38 - 1 : 46 ].strip()
    semimajorAxix = line[ 93 - 1 : 103 ].strip()
    orbitalEccentricity = line[ 71 - 1 : 79 ].strip()
    meanAnomalyEpoch = line[ 27 - 1 : 35 ].strip()
    slopeParamenter = line[ 15 - 1 : 19 ].strip()

    century = line[ 21 - 1 : 21 ].strip()
    lastTwoDigitsOfYear = line[ 22 - 1 : 23 ].strip()
    year = str( centuryMap[ century ] + int( lastTwoDigitsOfYear ) )
    month = getUnpackedDate( line[ 24 - 1 : 24 ].strip() )
    day = getUnpackedDate( line[ 25 - 1 : 25 ].strip() )
    epochDate = month + '/' + day + '/' + year

    xEphemLine = [
        name, 'e', inclinationToEcliptic, longitudeAscendingNode,
        argumentPerihelion, semimajorAxix, '0', orbitalEccentricity,
        meanAnomalyEpoch, epochDate, "2000.0", absoluteMagnitude, slopeParamenter ]

    xEphemLine = ','.join( xEphemLine )

    return xEphemLine


def getApparentMagnitudeOverPeriod( orbitalElement, startDate, endDate, city, apparentMagnitudeMaximum, intervalInDaysBetweenBodyCompute ):
    currentDate = startDate
    body = ephem.readdb( orbitalElement )
    while currentDate < endDate:
        [city.date](http://city.date/) = currentDate
        body.compute( city )
        if body.mag <= apparentMagnitudeMaximum:
            break

        currentDate += datetime.timedelta( days = intervalInDaysBetweenBodyCompute )

    return body.mag


# The parameters below take approximately 5 minutes to yield 1200 filtered bodies (as of 2022).
apparentMagnitudeMaximum = 15.0 # Upper limit for apparent magnitude.
durationInYears = 1 # Number of years over which to compute apparent magnitude.
intervalInDaysBetweenBodyCompute = 30

inFile = "MPCORB.DAT"
outFile = inFile[ 0 : -4 ] + "-FILTERED-BY-APPARENT-MAGNITUDE-" + str( int( apparentMagnitudeMaximum ) ) + ".DAT"
startDate = datetime.datetime.utcnow()
endDate = startDate + datetime.timedelta( days = 366 * durationInYears )
with open( inFile, 'r' ) as fIn, open( outFile, 'w' ) as fOut:
    city = ephem.city( "London" ) # Use any city; makes no difference to obtain the apparent magnitude.
    header = ""
    endOfHeader = False
    for line in fIn:
        if endOfHeader:
            if len( line.strip() ) == 0:
                continue

            name = line[ 167 - 1 : 194 ].replace( '(', '' ).replace( ')', '' ).strip()
            if len( name ) == 0:
                print( "Missing name:\n" + line )
                continue

            absoluteMagnitude = line[ 9 - 1 : 13 ].strip()
            if len( absoluteMagnitude ) == 0:
                print( "Missing absolute magnitude:\n" + line )
                continue

            apparentMagnitude = getApparentMagnitudeOverPeriod( mpcToXEphem( line ), startDate, endDate, city, apparentMagnitudeMaximum, intervalInDaysBetweenBodyCompute )
            if apparentMagnitude < apparentMagnitudeMaximum:
                fOut.write( line )

        else:
            header += line
            if line.startswith( "----------" ):
                endOfHeader = True
                fOut.write( header )

Conclusion

Unfortunately for my purposes, I have hit the wall. My astronomy program is an appindicator running on a user's computer and downloads the data files (bright/critical/distant/unusual) once a week. The files are small and so there is little overhead (to download and process).

It is now unreasonable to switch to download the MPCORB.DAT which is around 250MB once a week and then process. I did toy with the idea of downloading myself, generating the filtered version and then making available for download elsewhere, say GitHub, but I"m not sure how "proper" that would be.

Unless I can cross this (hopefully) last hurdle, I have to turn off the minor planet functionality. Note that this issue would still arise had I been fully loaded with Skyfield.

Thankfully the single data file for comets is small and matches up well with what I was originally using. See below for a port of the Perl script mentioned earlier which converts MPC format to XEphem/PyEphem:

Convert comets in MPC format to XEphem/PyEphem:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert comet data in MPC format to XEphem format.
#
# Inspired by:
#    https://github.com/XEphem/XEphem/blob/main/GUI/xephem/tools/mpccomet2edb.pl
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/CometOrbitFormat.html
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import sys


def processAndWriteOneLine( line, outputFile ):
    if len( line.strip() ) > 0:
        name = line[ 103 - 1 : 158 ].replace( '(', '' ).replace( ')', '' ).strip()
        absoluteMagnitude = line[ 92 - 1 : 95 ].strip() # $G The Perl script uses 91 instead of 92.

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            inclination = line[ 72 - 1 : 79 ].strip() # $i The Perl script uses 71 instead of 72. 
            longitudeAscendingNode = line[ 62 - 1 : 69 ].strip() # $O The Perl script uses 61 instead of 62.
            argumentPerihelion = line[ 52 - 1 : 59 ].strip() # $o The Perl script uses 51 instead of 52.
            perihelionDistance = line[ 31 - 1 : 39 ].strip() # $q
            orbitalEccentricity = line[ 42 - 1 : 49 ].strip() # $e The Perl script uses 41 instead of 42.

            month = line[ 20 - 1 : 21 ].strip()
            day = line[ 23 - 1 : 29 ].strip()
            year = line[ 15 - 1 : 18 ].strip()
            epochDate = month + '/' + day + '/' + year # $E

            slopeParameter = line[ 97 - 1 : 100 ].strip() # $H

            if float( orbitalEccentricity ) < 0.99: # Elliptical orbit.
                meanAnomaly = str( 0.0 ) # $M
                meanDistance = str( float( perihelionDistance ) / ( 1.0 - float( orbitalEccentricity ) ) ) # $a

                components = [
                    name, 'e', inclination, longitudeAscendingNode, argumentPerihelion,
                    meanDistance, '0', orbitalEccentricity, meanAnomaly,
                    epochDate, "2000.0",
                    absoluteMagnitude, slopeParameter ]

            elif float( orbitalEccentricity ) > 1.0: # Hyperbolic orbit.
                components = [
                    name, 'h', epochDate, inclination,
                    longitudeAscendingNode, argumentPerihelion, orbitalEccentricity, 
                    perihelionDistance, "2000.0",
                    absoluteMagnitude, slopeParameter ]

            else: # Parabolic orbit.
                components = [
                    name, 'p', epochDate, inclination,
                    argumentPerihelion, perihelionDistance, longitudeAscendingNode, "2000.0",
                    absoluteMagnitude, slopeParameter ]

            outputFile.write( ','.join( components ) + '\n' )


if len( sys.argv ) != 2:
    message = \
        "Usage: python3 mpcCometToXEphem.py fileToConvert" + \
        "\n\nFor example:" + \
        "\n  python3 mpcCometToXEphem.py CometEls.txt"

    raise SystemExit( message )

inFile = sys.argv[ 1 ]
outFile = inFile[ 0 : -3 ] + "edb"
with open( inFile, 'r' ) as fIn, open( outFile, 'w' ) as fOut:
    for line in fIn:
        processAndWriteOneLine( line, fOut )

[Bernmeister]

To clarify, the two scripts included above each serve a different purpose and had different constraints...

Convert comets in MPC format to XEphem/PyEphem: This script converts from MPC format to XEphem / PyEphem. No filtering by apparent magnitude and no dependency on PyEphem. This script can happily live in XEphem and/or PyEphem as a utility and essentially matches the Perl version. Given the comet data file size is quite small compared to MPCORB.DAT, I felt that filtering by apparent magnitude was not required and can easily be done by the end user within their own PyEphem code. I'm happy to add filtering (or create an additional script which does conversion and filtering).

Filter MPCORB.DAT by apparent magnitude: This script filters the MPCORB.DAT to a smaller file, leaving the smaller file in the same format as the original. The filtering is by apparent magnitude and uses PyEphem (rather than Skyfield due to #490). This script has a dependency on PyEphem (to annoy Skyfield users) and leaves the final filtered file in MPC format (to annoy PyEphem users). To my thinking, both Skyfield and PyEphem (and XEphem) users would benefit from a filtered version of MPCORB.DAT in their final respective formats. To that end, I'm happy to cobble together...

• A script which converts from MPCORB.DAT to XEphem / PyEphem. No filtering, no dependency on PyEphem and lives alongside the Perl equivalent. Somewhat useless on its own, but does just one job.

• A script which filters MPCORB.DAT and outputs to the same format (the filtering would be performed by PyEphem for now and Skyfield when up for it). Useful for Skyfield users.

• A script which filters MPCORB.DAT and outputs to the PyEphem / XEphem format (the filtering would be performed by PyEphem). Useful for PyEphem / XEphem users. Perhaps redundant given the script to filter and then the script to convert formats...

Finally, if feasible, running a script on GitHub to retrieve MPCORB.DAT, filter and provide for download, say weekly, would do the trick!

[Bernmeister]

Convert minor planets in MPC format to XEphem/PyEphem:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert minor planet data in MPC format to XEphem format.
#
# Inspired by:
#    https://github.com/XEphem/XEphem/blob/main/GUI/xephem/auxil/mpcorb2edb.pl
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import gzip, sys


centuryMap = {
    'I': 1800,
    'J': 1900,
    'K': 2000 }


# https://www.minorplanetcenter.net/iau/info/PackedDates.html
def getUnpackedDate( i ): return str( i ) if i.isdigit() else str( ord( i ) - ord( 'A' ) + 10 )


def processAndWriteOneLine( line, outputFile ):
    if len( line.strip() ) > 0:
        name = line[ 167 - 1 : 194 ].replace( '(', '' ).replace( ')', '' ).strip()
        absoluteMagnitude = line[ 9 - 1 : 13 ].strip() # $H

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            inclinationToEcliptic = line[ 60 - 1 : 68 ].strip() # $i
            longitudeAscendingNode = line[ 49 - 1 : 57 ].strip() # $O
            argumentPerihelion = line[ 38 - 1 : 46 ].strip() # $o
            semimajorAxix = line[ 93 - 1 : 103 ].strip() # $a
            orbitalEccentricity = line[ 71 - 1 : 79 ].strip() # $e
            meanAnomalyEpoch = line[ 27 - 1 : 35 ].strip() # $M
            slopeParamenter = line[ 15 - 1 : 19 ].strip() # $G

            century = line[ 21 - 1 : 21 ].strip() # $cent
            lastTwoDigitsOfYear = line[ 22 - 1 : 23 ].strip() # $TY
            year = str( centuryMap[ century ] + int( lastTwoDigitsOfYear ) ) # $TY
            month = getUnpackedDate( line[ 24 - 1 : 24 ].strip() ) # $TM
            day = getUnpackedDate( line[ 25 - 1 : 25 ].strip() ) # $TD
            epochDate = month + '/' + day + '/' + year

            components = [
                name, 'e', inclinationToEcliptic, longitudeAscendingNode,
                argumentPerihelion, semimajorAxix, '0', orbitalEccentricity, meanAnomalyEpoch,
                epochDate, "2000.0",
                absoluteMagnitude, slopeParamenter ]

            outputFile.write( ','.join( components ) + '\n' )


def convertTXT( inFile ):
    outFile = inFile[ 0 : -3 ] + "edb"
    with open( inFile, 'r' ) as fIn, open( outFile, 'w' ) as fOut:
        for line in fIn:
            processAndWriteOneLine( line, fOut )


def convertMPCORB( inFile ):
    outFile = inFile.replace( "DAT", "edb" )
    if inFile.endswith( ".gz" ):
        fIn = gzip.open( inFile, 'r' )
        fOut = gzip.open( outFile, 'w' )

    else:
        fIn = open( inFile, 'r' )
        fOut = open( outFile, 'w' )

    endOfHeader = False
    for line in fIn:
        if endOfHeader:
            processAndWriteOneLine( line, fOut )

        elif line.startswith( "----------" ):
            endOfHeader = True

    fIn.close()
    fOut.close()


if len( sys.argv ) != 2:
    message = \
        "Usage: python3 mpcMinorPlanetToXEphem.py fileToConvert" + \
        "\n\nFor example:" + \
        "\n  python3 mpcMinorPlanetToXEphem.py MPCORB.DAT" + \
        "\n  python3 mpcMinorPlanetToXEphem.py MPCORB.DAT.gz" + \
        "\n  python3 mpcMinorPlanetToXEphem.py mpcDataAsStraightTextWithoutHeader.txt"

    raise SystemExit( message )

inFile = sys.argv[ 1 ]

if inFile.endswith( ".txt" ):
    convertTXT( inFile )

elif inFile.startswith( "MPCORB.DAT" ):
    convertMPCORB( inFile )

else:
    raise SystemExit( "Unknown input file format." )

[Bernmeister]

I have been calculating the apparent magnitude of comets / minor planets using both PyEphem and Skyfield, comparing the results with each other and against some external sources. I have come to the conclusion that the magnitude parameters in the MPC data files are somewhat dated. Refer to #416 and the comment I made on 17 June 2022.

For those following along, I have found external sources which give (what I believe) is a more accurate/current value of apparent magnitude. Therefore, the whole issue around filtering MPCORB.DAT and even CometEls.txt is somewhat moot at this point.

More to come...

[rmathar]

The perl script mentioned above, GUI/xephem/auxil/mpcorb2edb.pl, takes the entire MPCORB.DAT file as input and generates two files, filtering the objects by magnitude. The threshold of the limiting magnitude is the parameter dimmag in (roughly) line 263 of the perl code. Because the apparent magnitude of the asteroids is obviously a function of time (depending on the distance sun-to-planet and planet-to-earth etc), the limiting magnitude is roughly estimated by the periapsis/apoapsis of the planet's orbit.

[Bernmeister]

The perl script mentioned above, GUI/xephem/auxil/mpcorb2edb.pl, takes the entire MPCORB.DAT file as input and generates two files, filtering the objects by magnitude. The threshold of the limiting magnitude is the parameter dimmag in (roughly) line 263 of the perl code. Because the apparent magnitude of the asteroids is obviously a function of time (depending on the distance sun-to-planet and planet-to-earth etc), the limiting magnitude is roughly estimated by the periapsis/apoapsis of the planet's orbit.

The problem (as I understand it) is (for comets at least) the absolute magnitude is out of date and using that to then compute the apparent magnitude is erroneous.

[Bernmeister]

Epilogue: I have found a source for each of minor planets and comets which provide updated orbital elements and corresponding apparent magnitude (either within the orbital elements or as a separate piece of data).

For minor planets, I am using data from Lowell Minor Planet Services; for comets, the Comet Observation database.

In terms of the scripts above which convert/filter, please ignore those, as I've found errors and changed my mind about use cases. Instead, I have put together replacement Python3 scripts which implement:

• Lowell minor planet to MPC
• Lowell minor planet to XEphem
• MPC comet to XEphem
• MPC minorplanet to XEphem

These scripts are in the same vein as the Perl scripts in the XEphem repository mentioned above.

@brandon-rhodes I don't want to edit the previous posts/replies to remove old scripts in an attempt to "tidy up". I am happy to paste in here my new scripts and they are free for all. If you want to include the scripts in Skyfield, PyEphem or XEphem repositories, you have my blessing.

[brandon-rhodes]

@Bernmeister — I would be happy to see your new scripts, if you would paste them in. Thanks!

[Bernmeister]

@brandon-rhodes Happy for these to be incorporated into the PyEphem (or XEphem) repository, at your discretion.

Convert minor planet data in Lowell format to MPC format

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert minor planet data in Lowell format to MPC format.
#
# Lowell format: 
#    https://asteroid.lowell.edu/main/astorb/
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html


import argparse, gzip


# https://www.minorplanetcenter.net/iau/info/PackedDates.html
def getPackedDate( year, month, day ):
    packedYear = year[ 2 : ]
    if int( year ) < 1900:
        packedYear = 'I' + packedYear

    elif int( year ) < 2000:
        packedYear = 'J' + packedYear

    else:
        packedYear = 'K' + packedYear


    def getPackedDayMonth( dayOrMonth ):
        if int( dayOrMonth ) < 10:
            packedDayMonth = str( int( dayOrMonth ) )

        else:
            packedDayMonth = chr( int( dayOrMonth ) - 10 + ord( 'A' ) )

        return packedDayMonth

    packedMonth = getPackedDayMonth( month )
    packedDay = getPackedDayMonth( day )

    return packedYear + packedMonth + packedDay


def processAndWriteOneLine( line, outputFile ):
    if len( line.strip() ) > 0:
        startIndices = [ 1, 8, 27, 43, 49, 55, 60, 66, 71, 96, 101, 107, 116, 127, 138, 148, 159, 169, 183, 192, 200, 209, 218, 235, 252 ]
        startIndices = [ x - 1 for x in startIndices ] # Offset back to zero to match each line read into a string.

        parts = [ "OFFSET TO ALIGN WITH FIELD NUMBERING" ] + \
                [ line[ i : j ] for i, j in zip( startIndices, startIndices[ 1 : ] + [ None ] ) ] # Inspired by https://stackoverflow.com/a/10851479/2156453

        number = parts[ 1 ].strip()
        name = parts[ 2 ].strip()
        absoluteMagnitude = parts[ 4 ].strip()

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            slopeParameter = parts[ 5 ].strip()
            numberObservations = parts[ 11 ].strip()
            epochDate = parts[ 12 ].strip()
            meanAnomalyEpoch = parts[ 13 ].strip()
            argumentPerihelion = parts[ 14 ].strip()
            longitudeAscendingNode = parts[ 15 ].strip()
            inclinationToEcliptic = parts[ 16 ].strip()
            orbitalEccentricity = parts[ 17 ].strip()
            semimajorAxis = parts[ 18 ].strip()

            components = [
                ' ' * 7, # number or designation packed
                ' ', # 8
                str( round( float( absoluteMagnitude ), 2 ) ).rjust( 5 ),
                ' ', # 14
                str( round( float( slopeParameter ), 2 ) ).rjust( 5 ),
                ' ', # 20
                getPackedDate( epochDate[ 0 : 4 ], epochDate[ 4 : 6 ], epochDate[ 6 : 8 ] ).rjust( 5 ),
                ' ', # 26
                str( round( float( meanAnomalyEpoch ), 5 ) ).rjust( 9 ),
                ' ' * 2, # 36, 37 
                str( round( float( argumentPerihelion ), 5 ) ).rjust( 9 ),
                ' ' * 2, # 47, 48
                str( round( float( longitudeAscendingNode ), 5 ) ).rjust( 9 ),
                ' ' * 2, # 58, 59
                str( round( float( inclinationToEcliptic ), 5 ) ).rjust( 9 ),
                ' ' * 2, # 69, 70
                str( round( float( orbitalEccentricity ), 7 ) ).rjust( 9 ),
                ' ', # 80
                ' ' * 11, # mean daily motion
                ' ', # 92
                str( round( float( semimajorAxis ), 7 ) ).rjust( 11 ),
                ' ' * 2, # 104, 105
                ' ', # uncertainty parameter
                ' ', # 107
                ' ' * 9, # reference
                ' ', # 117
                numberObservations.rjust( 5 ),
                ' ', # 123
                ' ' * 3, # oppositions
                ' ', # 127
                ' ' * ( 4 + 1 + 4 ), # multiple/single oppositions
                ' ', # 137
                ' ' * 4, # rms residual
                ' ', # 142
                ' ' * 3, # coarse indicator of perturbers
                ' ', # 146
                ' ' * 3, # precise indicator of perturbers
                ' ', # 150
                ' ' * 10, # computer name
                ' ', # 161
                ' ' * 4, # hexdigit flags
                ' ', # 166
                ( number + ' ' + name ).strip().ljust( 194 - 167 + 1 ),
                ' ' * 8 ] # date last observation

            outputFile.write( ''.join( components ) + '\n' )


def convert( inFile, outFile ):
    if inFile.endswith( ".gz" ):
        fIn = gzip.open( inFile, 'rt' )

    else:
        fIn = open( inFile, 'r' )

    fOut = open( outFile, 'w' )
    for line in fIn:
        processAndWriteOneLine( line, fOut )

    fIn.close()
    fOut.close()


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description = "Convert a minor planet text file such as astorb.dat or astorb.dat.gz from Lowell to MPC format. " + \
                      "If the file ends in '.gz' the file will be treated as a gzip file; otherwise the file is assumed to be text." )
    parser.add_argument( "inFile", help = "File to convert" )
    parser.add_argument( "outFile", help = "Output file to be created" )
    args = parser.parse_args()
    convert( args.inFile, args.outFile )

Convert minor planet data in Lowell format to XEphem format

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert minor planet data in Lowell format to XEphem format.
#
# Inspired by:
#    https://github.com/XEphem/XEphem/blob/main/GUI/xephem/auxil/astorb2edb.pl
#
# Lowell format: 
#    https://asteroid.lowell.edu/main/astorb/
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import argparse, gzip


def processAndWriteOneLine( line, outputFile ):
    if len( line.strip() ) > 0:
        startIndices = [ 1, 8, 27, 43, 49, 55, 60, 66, 71, 96, 101, 107, 116, 127, 138, 148, 159, 169, 183, 192, 200, 209, 218, 235, 252 ]
        startIndices = [ x - 1 for x in startIndices ] # Offset back to zero to match each line read into a string.

        parts = [ "OFFSET TO ALIGN WITH FIELD NUMBERING" ] + \
                [ line[ i : j ] for i, j in zip( startIndices, startIndices[ 1 : ] + [ None ] ) ] # Inspired by https://stackoverflow.com/a/10851479/2156453

        number = parts[ 1 ].strip()
        name = parts[ 2 ].strip()
        absoluteMagnitude = parts[ 4 ].strip()

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            slopeParameter = parts[ 5 ].strip()
            epochDate = parts[ 12 ].strip()
            meanAnomalyEpoch = parts[ 13 ].strip()
            argumentPerihelion = parts[ 14 ].strip()
            longitudeAscendingNode = parts[ 15 ].strip()
            inclinationToEcliptic = parts[ 16 ].strip()
            orbitalEccentricity = parts[ 17 ].strip()
            semimajorAxis = parts[ 18 ].strip()

            components = [
                ( number + ' ' + name ).strip(),
                'e',
                inclinationToEcliptic,
                longitudeAscendingNode,
                argumentPerihelion,
                semimajorAxis,
                '0',
                orbitalEccentricity,
                meanAnomalyEpoch,
                epochDate[ 4 : 6 ] + '/' + epochDate[ 6 : ] + '/' + epochDate[ 0 : 4 ],
                "2000.0",
                absoluteMagnitude,
                slopeParameter ]

            outputFile.write( ','.join( components ) + '\n' )


def convert( inFile, outFile ):
    if inFile.endswith( ".gz" ):
        fIn = gzip.open( inFile, 'rt' )

    else:
        fIn = open( inFile, 'r' )

    fOut = open( outFile, 'w' )
    for line in fIn:
        processAndWriteOneLine( line, fOut )

    fIn.close()
    fOut.close()


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description = "Convert a minor planet text file such as astorb.dat or astorb.dat.gz from Lowell to XEphem format. " + \
                      "If the file ends in '.gz' the file will be treated as a gzip file; otherwise the file is assumed to be text." )
    parser.add_argument( "inFile", help = "File to convert" )
    parser.add_argument( "outFile", help = "Output file to be created" )
    args = parser.parse_args()
    convert( args.inFile, args.outFile )

Convert comet data in MPC format to XEphem format

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert comet data in MPC format to XEphem format.
#
# Inspired by:
#    https://github.com/XEphem/XEphem/blob/main/GUI/xephem/tools/mpccomet2edb.pl
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/CometOrbitFormat.html
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import argparse


def processAndWriteOneLine( line, outFile ):
    if len( line.strip() ) > 0:
        # Field numbers: 0  1   2   3   4   5   6   7   8   9  10  11  12  13  14   15   16   17
        startIndices = [ 1, 5,  6, 15, 20, 23, 31, 42, 52, 62, 72, 82, 86, 88, 92,  97, 103, 160 ]
        endIndices =   [ 4, 5, 12, 18, 21, 29, 39, 49, 59, 69, 79, 85, 87, 89, 95, 100, 158, 168 ]

        # Offset back to zero to match lines read into a string.
        startIndices = [ x - 1 for x in startIndices ]
        endIndices = [ x - 1 for x in endIndices ]

        fields = [ line[ i : j + 1 ] for i, j in zip( startIndices, endIndices ) ] # Inspired by https://stackoverflow.com/a/10851479/2156453
        name = fields[ 16 ].replace( '(', '' ).replace( ')', '' ).strip()
        absoluteMagnitude = fields[ 14 ].strip()

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            year = fields[ 3 ].strip()
            month = fields[ 4 ].strip()
            day = fields[ 5 ].strip()
            epochDate = month + '/' + day + '/' + year

            perihelionDistance = fields[ 6 ].strip()
            orbitalEccentricity = fields[ 7 ].strip()
            argumentPerihelion = fields[ 8 ].strip()
            longitudeAscendingNode = fields[ 9 ].strip()
            inclination = fields[ 10 ].strip()
            slopeParameter = fields[ 15 ].strip()

            if float( orbitalEccentricity ) < 0.99: # Elliptical orbit.
                meanAnomaly = str( 0.0 )
                meanDistance = str( float( perihelionDistance ) / ( 1.0 - float( orbitalEccentricity ) ) )

                components = [
                    name,
                    'e',
                    inclination,
                    longitudeAscendingNode,
                    argumentPerihelion,
                    meanDistance,
                    '0',
                    orbitalEccentricity,
                    meanAnomaly,
                    epochDate,
                    "2000.0",
                    absoluteMagnitude,
                    slopeParameter ]

            elif float( orbitalEccentricity ) > 1.0: # Hyperbolic orbit.
                components = [
                    name,
                    'h',
                    epochDate,
                    inclination,
                    longitudeAscendingNode,
                    argumentPerihelion,
                    orbitalEccentricity, 
                    perihelionDistance,
                    "2000.0",
                    absoluteMagnitude,
                    slopeParameter ]

            else: # Parabolic orbit.
                components = [
                    name,
                    'p',
                    epochDate,
                    inclination,
                    argumentPerihelion,
                    perihelionDistance,
                    longitudeAscendingNode,
                    "2000.0",
                    absoluteMagnitude,
                    slopeParameter ]

            outFile.write( ','.join( components ) + '\n' )


def convert( inFile, outFile ):
    fIn = open( inFile, 'r' )
    fOut = open( outFile, 'w' )
    for line in fIn:
        processAndWriteOneLine( line, fOut )

    fIn.close()
    fOut.close()


if __name__ == "__main__":
    parser = argparse.ArgumentParser( description = "Convert a comet text file such as CometEls.txt from MPC to XEphem format." )
    parser.add_argument( "inFile", help = "File to convert" )
    parser.add_argument( "outFile", help = "Output file to be created" )
    args = parser.parse_args()
    convert( args.inFile, args.outFile )

Convert minor planet data in MPC format to XEphem format

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.


# Convert minor planet data in MPC format to XEphem format.
#
# Inspired by:
#    https://github.com/XEphem/XEphem/blob/main/GUI/xephem/auxil/mpcorb2edb.pl
#
# MPC format: 
#    https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html
#
# XEphem format:
#    https://xephem.github.io/XEphem/Site/help/xephem.html#mozTocId468501


import argparse, gzip


centuryMap = {
    'I': 1800,
    'J': 1900,
    'K': 2000 }


# https://www.minorplanetcenter.net/iau/info/PackedDates.html
def getUnpackedDate( i ): return str( i ) if i.isdigit() else str( ord( i ) - ord( 'A' ) + 10 )


def processAndWriteOneLine( line, outputFile ):
    if len( line.strip() ) > 0:
        # Field numbers: 0   1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17   18   19   20   21   22
        startIndices = [ 1,  9, 15, 21, 27, 38, 49, 60, 71, 81,  93, 106, 108, 118, 124, 128, 138, 143, 147, 151, 162, 167, 195 ]
        endIndices =   [ 7, 13, 19, 25, 35, 46, 57, 68, 79, 91, 103, 106, 116, 122, 126, 136, 141, 145, 149, 160, 165, 194, 202 ]

        # Offset back to zero to match lines read into a string.
        startIndices = [ x - 1 for x in startIndices ]
        endIndices = [ x - 1 for x in endIndices ]

        fields = [ line[ i : j + 1 ] for i, j in zip( startIndices, endIndices ) ] # Inspired by https://stackoverflow.com/a/10851479/2156453
        name = fields[ 21 ].replace( '(', '' ).replace( ')', '' ).strip()
        absoluteMagnitude = fields[ 1 ].strip()

        if len( name ) == 0:
            print( "Missing name:\n" + line )

        elif len( absoluteMagnitude ) == 0:
            print( "Missing absolute magnitude:\n" + line )

        else:
            slopeParameter = fields[ 2 ].strip()

            epochPacked = fields[ 3 ].strip()
            century = epochPacked[ 0 ]
            lastTwoDigitsOfYear = epochPacked[ 1 : 3 ]
            year = str( centuryMap[ century ] + int( lastTwoDigitsOfYear ) )
            month = getUnpackedDate( epochPacked[ 3 ] )
            day = getUnpackedDate( epochPacked[ 4 ] )
            epochDate = month + '/' + day + '/' + year

            meanAnomalyEpoch = fields[ 4 ].strip()
            argumentPerihelion = fields[ 5 ].strip()
            longitudeAscendingNode = fields[ 6 ].strip()
            inclinationToEcliptic = fields[ 7 ].strip()
            orbitalEccentricity = fields[ 8 ].strip()
            semimajorAxis = fields[ 10 ].strip()

            components = [
                name,
                'e',
                inclinationToEcliptic,
                longitudeAscendingNode,
                argumentPerihelion,
                semimajorAxis,
                '0',
                orbitalEccentricity,
                meanAnomalyEpoch,
                epochDate,
                "2000.0",
                absoluteMagnitude,
                slopeParameter ]

            outputFile.write( ','.join( components ) + '\n' )


def convert( inFile, header, outFile ):
    if inFile.endswith( ".gz" ):
        fIn = gzip.open( inFile, 'rt' )

    else:
        fIn = open( inFile, 'r' )

    fOut = open( outFile, 'w' )

    if header:
        endOfHeader = False
        for line in fIn:
            if endOfHeader:
                processAndWriteOneLine( line, fOut )
    
            elif line.startswith( "----------" ):
                endOfHeader = True

    else:
        for line in fIn:
            processAndWriteOneLine( line, fOut )

    fIn.close()
    fOut.close()


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description = "Convert a minor planet text file such as MPCORB.DAT or MPCORB.DAT.gz from MPC to XEphem format. " + \
        "If the file ends in '.gz' the file will be treated as a gzip file; otherwise the file is assumed to be text." )
    parser.add_argument( "inFile", help = "File to convert" )
    parser.add_argument(
        "--header",
        action = "store_true",
        help = "Only specify this option if the file to convert contains a header, such as MPCORB.DAT or MPCORB.DAT.gz." )
    parser.add_argument( "outFile", help = "Output file to be created" )
    args = parser.parse_args()
    convert( args.inFile, args.header, args.outFile )