comets and mps

[heinkurz]

hi! how can i get the current 10 brightest comets and minor planets?

[EndlessDex]

Here is an example that should get you close.

https://rhodesmill.org/skyfield/stars.html#filtering-the-star-catalog

Note that what results you get will greatly depend on which catalog you use.

[heinkurz]

sorry but how could stars help to find comets and minor planets ???

[EndlessDex]

Sorry for not being clearer. You use the same filtering process with a different catalog. Probably the "mpc" catalog.

[EndlessDex]

Here is an example of fetching data from MPC

https://rhodesmill.org/skyfield/example-plots.html#drawing-a-finder-chart-for-comet-neowise

[albertw]

You probably want to grab https://minorplanetcenter.net/iau/Ephemerides/Bright/2025/Soft00Bright.txt

Then have a look through #416 which seems to discus methods of calculating the apparent magnitude

[EndlessDex]

Apparent vs absolute 😵‍💫

[heinkurz]

I'm now near a solution and it works with comets but I have a problem with minor planets.

    cdf = pd.DataFrame()
    with api.load.open("MPCORB10k.DAT") as f:       # reduced file of MPCORB.DAT with the first 10k minor planets
        mpcs = mpc.load_mpcorb_dataframe(f)
    mpcs = (mpcs.sort_values('reference').groupby('designation', as_index=False).last().set_index('designation', drop=False))
    ts = api.load.timescale()
    eph = api.load('de421.bsp')
    sun = eph['sun']

    for index, row in mpcs.iterrows():
        mp = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)
        .........

The last line raises an error:
p = a * (1.0 - e*e)
^
TypeError: can't multiply sequence by non-int of type 'str'

What's wrong? I've downloaded MPCORB.DAT from the MPC site so it must be correct.

[albertw]

Just a guess but perhaps this:
mpcs = (mpcs.sort_values('reference').groupby('designation', as_index=False).last().set_index('designation', drop=False))

is somehow setting values to strings that should be floats.

[heinkurz]

Hm, when I delete that line I get a different error in the same line:
year = 100 * n(s[0]) + int(s[1:3])
^^^^^^^^^^^
ValueError: invalid literal for int() with base 10: 'RB'

I think that means the data file format is wrong - or mpc.mpcorb_orbit() has a problem !?!?!?

[heinkurz]

Problem found!
MPCORB.DAT must be free of all header lines, only the pure data lines should remain.
And there have been failures in the original source code.
So here is my solution:

import datetime
import math
import pytz
import pandas as pd
from skyfield import api
from skyfield.data import mpc
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN

local = pytz.timezone("UTC")
today = datetime.date.today()
start_date = datetime.datetime(today.year, today.month, today.day)
start_date = local.localize(start_date)
end_date = datetime.datetime(today.year, today.month, today.day + 1)
end_date = local.localize(end_date)

location = {
    "location": "Vienna",
    "location_lat": 48.20849,
    "location_long": 16.37208,
    "location_elev": 200.0,
    "tz": pytz.timezone("Europe/Vienna")
}

lim_mag = 10.0
lim_sun_elongation = 10
lim_moon_elongation = 10



# ######### comets #########

def comets():

    cdf = pd.DataFrame()

    with api.load.open(mpc.COMET_URL) as f:
        comets = mpc.load_comets_dataframe_slow(f)

    comets = (comets.sort_values("reference").groupby("designation", as_index=False).last().set_index("designation", drop=False))

    ts = api.load.timescale()
    eph = api.load("de421.bsp")
    sun = eph["sun"]
    moon = eph["moon"]
    earth = eph["earth"]
    loc = api.wgs84.latlon(location["location_lat"], location["location_long"])
    observer = earth + api.Topos(location["location_lat"], location["location_long"],
                                          elevation_m=location["location_elev"])
    t0 = ts.utc(start_date)
    observed_sun = observer.at(t0).observe(sun)
    observed_moon = observer.at(t0).observe(moon)

    z = 0
    for index, row in comets.iterrows():
        comet = sun + mpc.comet_orbit(row, ts, GM_SUN)
        observed_comet = observer.at(t0).observe(comet)

        r = comet.at(t0).distance()
        ra, dec, delta = observed_comet.radec()
        sun_elong = observed_sun.separation_from(observed_comet).degrees
        moon_elong = observed_moon.separation_from(observed_comet).degrees
        mag = row.magnitude_g + 5 * math.log10(delta.au) + 2.5 * row.magnitude_k * math.log10(r.au)

        c = pd.DataFrame({
             "name": row.designation,
             "ra": ra.hours,
             "dec": dec.degrees,
             "delta": delta.au,
             "sun_elongation": sun_elong,
             "moon_elongation": moon_elong,
             "magnitude": mag}, index=[0]
            )
        #cdf = pd.concat([c, cdf])
        if mag <= lim_mag:
            print(c.to_string() + "\r\n")
            z += 1

    #cdf = cdf.loc[(cdf["magnitude"] <= lim_mag) & (cdf["sun_elongation"] >= lim_sun_elongation) & (cdf["moon_elongation"] >= lim_moon_elongation)]
    #print(cdf.to_string())
    print(z, "found")



# ######### minor planets #########

def minorplanets():

    cdf = pd.DataFrame()

    #with api.load.open(mpc.MPCORB_URL) as f:   # , decoding="utf-8"
    with api.load.open("MPCORB10k.DAT") as f:
        mpcs = mpc.load_mpcorb_dataframe(f)

    #mpcs = (mpcs.sort_values("reference").groupby("designation",
    #                                              as_index=False).last().set_index("designation", drop=False))
    #print(mpcs)

    ts = api.load.timescale()
    eph = api.load("de421.bsp")
    sun = eph["sun"]
    moon = eph["moon"]
    earth = eph["earth"]
    loc = api.wgs84.latlon(location["location_lat"], location["location_long"])
    observer = earth + api.Topos(location["location_lat"], location["location_long"],
                                          elevation_m=location["location_elev"])
    t0 = ts.utc(start_date)
    observed_sun = observer.at(t0).observe(sun)
    observed_moon = observer.at(t0).observe(moon)

    z = 0
    for index, row in mpcs.iterrows():
        mp = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)
        observed_mp = observer.at(t0).observe(mp)
        r = mp.at(t0).distance()
        ra, dec, delta = observed_mp.radec()
        sun_elong = observed_sun.separation_from(observed_mp).degrees
        moon_elong = observed_moon.separation_from(observed_mp).degrees
        mag = row.magnitude_H + 5 * math.log10(delta.au) + 2.5 * row.magnitude_G * math.log10(r.au)

        c = pd.DataFrame({
             "name": row.designation,
             "ra": ra.hours,
             "dec": dec.degrees,
             "delta": delta.au,
             "sun_elongation": sun_elong,
             "moon_elongation": moon_elong,
             "magnitude": mag}, index=[0]
            )
        #cdf = pd.concat([c, cdf])
        if mag <= lim_mag:
            print(c.to_string() + "\r\n")
            z += 1

    #cdf = cdf.loc[(cdf["magnitude"] <= lim_mag) & (cdf["sun_elongation"] >= lim_sun_elongation) & (cdf["moon_elongation"] >= lim_moon_elongation)]
    #print(cdf.to_string())
    print(z, "found")


comets()
minorplanets()