Currently visible comets from a location (observe multiple comets)

[jurezakrajsek]

Im preparing a comet planner for my website, where i would need to calulate the visibility of ~2000 comets from MPC list (or my database list) and eliminate them by brightnes, altitude above horizon (at transit time), apparent distance from the sun and moon , and calculate the time the comets is highest in the sky ...
It should work in real-time as the users are supposed to be able to change the date that the list should be produced for.

Currently im doing this with pyephem, which seems to work quite vell but for some comets issues a near-parrabolic error.

So the idea was to swith to skyfields
I can loop through the objects and observe() each individaully, but this process takes almost 60s for 2000 comets, and I didnt even start to calculate the transit time, or time best visible.

Here is part of the code i came up with so far:

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

def main():
    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': 'Crni Vrh Observatory',
        'location_lat': 45.94583,
        'location_long': 14.07111,
        'location_elev': 726.0,
        'tz': pytz.timezone('Europe/Paris')
    }

    lim_mag = 15
    lim_sun_elongation = 10
    lim_moon_elongation = 10

    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))

    print(comets)

    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)

    for index, row in comets.iterrows():
        print(row.designation)
        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])

    cdf = cdf.loc[(cdf['magnitude'] <= lim_mag) & (cdf['sun_elongation'] >= lim_sun_elongation) & (cdf['moon_elongation'] >= lim_moon_elongation)]
    print(cdf)


if __name__ == "__main__":
    main()

Is there any better way on how to do this calculation with skyfields on multiple comets at once?
Any ideas or directions would be appreciated.

Best regards
Jure