[astorb2ephem] Convert the astorb.dat database into an ephemeris database for a given epoch #orbital_elements #astorb #pyephem #ephemeris

astorb2ephem.py

#!/usr/local/bin/python
# encoding: utf-8
"""
*Convert the astorb.dat database into an ephemeris database for a given epoch (a solar-system snapshot). Note this script uses PyEphem to generate the ephemerides, so although it is fast, the resulting ephemerides are far from accuracte*

:Author:
    David Young

:Date Created:
    October  6, 2017

Usage:
    astorb2ephem <pathToAstorbDB> <pathToEphemDB> <mjd> <magLimit>

    pathToAstorbDB        path to the astorb input database file
    pathToEphemDB         path to save the resulting ephemeris database to
    mjd                   the epoch at which to generate the ephemeris database
    magLimit              a lower mag limit above which to report ephemerides

Options:
    -h, --help            show this help message
    -v, --version         show version
    -s, --settings        the settings file
"""
################# GLOBAL IMPORTS ####################
import sys
import os
from fundamentals import tools
import ephem
import codecs
import math
import re
import numpy as np
import healpy as hp
from fundamentals.mysql import convert_dictionary_to_mysql_table
from fundamentals.mysql import readquery
from fundamentals.renderer import list_of_dictionaries
from fundamentals.mysql import directory_script_runner
import datetime
import codecs


def main(arguments=None):
    """
    *The main function used when ``astorb2xephem.py`` is run as a single script from the cl*
    """

    # SETUP THE COMMAND-LINE UTIL SETTINGS
    su = tools(
        arguments=arguments,
        docString=__doc__,
        logLevel="WARNING",
        options_first=False,
        projectName=False
    )
    arguments, settings, log, dbConn = su.setup()

    # UNPACK REMAINING CL ARGUMENTS USING `EXEC` TO SETUP THE VARIABLE NAMES
    # AUTOMATICALLY
    for arg, val in arguments.iteritems():
        if arg[0] == "-":
            varname = arg.replace("-", "") + "Flag"
        else:
            varname = arg.replace("<", "").replace(">", "")
        if isinstance(val, str) or isinstance(val, unicode):
            exec(varname + " = '%s'" % (val,))
        else:
            exec(varname + " = %s" % (val,))
        if arg == "--dbConn":
            dbConn = val
        log.debug('%s = %s' % (varname, val,))

    # CONSTANTS
    magLimit = float(magLimit)
    nside = 1024
    pi = (4 * math.atan(1.0))
    DEG_TO_RAD_FACTOR = pi / 180.0
    RAD_TO_DEG_FACTOR = 180.0 / pi

    # OPEN AND READ THE astorb.dat FILE
    pathToReadFile = pathToAstorbDB
    try:
        log.debug("attempting to open the file %s" % (pathToReadFile,))
        readFile = codecs.open(pathToReadFile, encoding='utf-8', mode='r')
        thisData = readFile.read()
        readFile.close()
    except IOError, e:
        message = 'could not open the file %s' % (pathToReadFile,)
        log.critical(message)
        raise IOError(message)
    readFile.close()

    # THE PYEPHEM OBSERVER
    obs = ephem.Observer()
    obs.date = float(mjd) - 15019.5
    objects = []

    lines = thisData.split("\n")
    for l in lines:

        if len(l) < 50:
            continue
        d = {}
        ephemDict = {}

        # PARSE THE LINE FROM astorb.dat (UNNEEDED VALUES COMMENTED OUT)
        d["mpc_number"] = l[0:7].strip()
        d["name"] = l[7:26].strip()
        d["discoverer"] = l[26:41].strip()
        d["H_abs_mag"] = l[41:48].strip()
        d["G_slope"] = l[48:54].strip()
        d["color_b_v"] = l[54:59].strip()
        d["diameter_km"] = l[59:65].strip()
        d["class"] = l[65:71].strip()
        # d["int1"] = l[71:75].strip()
        # d["int2"] = l[75:79].strip()
        # d["int3"] = l[79:83].strip()
        # d["int4"] = l[83:87].strip()
        # d["int5"] = l[87:91].strip()
        # d["int6"] = l[91:95].strip()
        d["orbital_arc_days"] = l[95:101].strip()
        d["number_obs"] = l[101:106].strip()
        d["epoch"] = l[106:115].strip()
        d["M_mean_anomaly_deg"] = l[115:126].strip()
        d["o_arg_peri_deg"] = l[126:137].strip()
        d["O_long_asc_node_deg"] = l[137:148].strip()
        d["i_inclination_deg"] = l[148:158].strip()
        d["e_eccentricity"] = l[158:169].strip()
        d["a_semimajor_axis"] = l[169:182].strip()
        d["orbit_comp_date"] = l[182:191].strip()
        d["ephem_uncertainty_arcsec"] = l[191:199].strip()
        d["ephem_uncertainty_change_arcsec_day"] = l[199:208].strip()
        d["ephem_uncertainty_date"] = l[208:217].strip()
        # d["peak_ephem_uncertainty_next_arcsec"] = l[217:225].strip()
        # d["peak_ephem_uncertainty_next_date"] = l[225:234].strip()
        # d["peak_ephem_uncertainty_10_yrs_from_ceu_arcsec"] = l[
        #     217:225].strip()
        # d["peak_ephem_uncertainty_10_yrs_from_ceu_date"] = l[
        #     242:251].strip()
        # d["peak_ephem_uncertainty_10_yrs_from_peu_arcsec"] = l[
        #     251:259].strip()
        # d["peak_ephem_uncertainty_10_yrs_from_peu_date"] = l[
        #     259:].strip()

        yyyy = int(d["epoch"][:4])
        mm = int(d["epoch"][4:6])
        dd = int(d["epoch"][6:])
        d["epoch_xeph"] = "%(mm)s/%(dd)s/%(yyyy)s" % locals()
        xephemStr = "%(mpc_number)s %(name)s,e,%(i_inclination_deg)s,%(O_long_asc_node_deg)s,%(o_arg_peri_deg)s,%(a_semimajor_axis)s,0,%(e_eccentricity)s,%(M_mean_anomaly_deg)s,%(epoch_xeph)s,2000.0,%(H_abs_mag)s,%(G_slope)s" % d
        xephemStr = xephemStr.strip()

        # TIDYUP
        if len(d["mpc_number"]) == 0:
            d["mpc_number"] = None

        # PyEphem works in DUBLIN JD, to convert from MJD subtract 15019.5
        minorPlanet = ephem.readdb(xephemStr)
        minorPlanet.compute(obs)

        # CONVERT TO A TOPCENTRIC POSITION

        # http://rhodesmill.org/pyephem/radec
        # a_ra, a_dec - Astrometric Geocentric Position for the star atlas epoch you've specified
        # g_ra, g_dec - Apparent Geocentric Position for the epoch-of-date
        # ra, dec - Apparent Topocentric Position for the epoch-of-date
        if minorPlanet.mag > magLimit:
            continue

        ephemDict["ra_deg"] = minorPlanet.a_ra * RAD_TO_DEG_FACTOR
        ephemDict["dec_deg"] = minorPlanet.a_dec * RAD_TO_DEG_FACTOR
        ephemDict["mpc_number"] = d["mpc_number"]
        ephemDict["object_name"] = d["name"]
        try:
            ephemDict["healpix"] = hp.ang2pix(
                nside, theta=ephemDict["ra_deg"], phi=ephemDict["dec_deg"], lonlat=True)
        except:
            continue

        objects.append(ephemDict)

    filepath = pathToEphemDB
    dataSet = list_of_dictionaries(
        log=log,
        listOfDictionaries=objects
    )

    csv = dataSet.csv(filepath=pathToEphemDB)

    return

if __name__ == '__main__':
    main()