In [1]:
!which python; python -V;
from astropy import units as u

from perylune.orbit_tools import *
from perylune.interplanetary import *
from perylune.mpc import *

import time

/home/thomson/devel/perylune/venv/bin/python
Python 3.8.5
In [2]:
# First step is to load the orbit data. parse_txt() will load the orbital elements from text file in MPC format.
# This will be returned as an array. Each entry in that array has a list of 8 parameters or so.
elements = parse_txt("../data/mpc/mpcorb_extended.dat", limit=0, skip="------------------")

print("Number of orbital elements parsed: %d" % len(elements))

Read 1028498 line(s) from NEA.txt, parsing up to 0 line(s)
WARNING: Skipped line 546847: too short (1, expected at least 167 chars)
WARNING: Skipped line 978222: too short (1, expected at least 167 chars)
Number of orbital elements parsed: 1028453
In [3]:
# Now process the data. Find min and max a and e, also put the a and e values in their separate arrays (a_tbl, e_tbl)

max_e = 0
min_e = 9999999
max_a = 0
min_a = 9999999

a_tbl = []
e_tbl = []
cnt = 0
for e in elements:
    a = e[8].value
    ecc = e[6].value
    if a < 6:
        cnt = cnt + 1
        a_tbl.append(a)
        e_tbl.append(e)
    if max_a < a:
        max_a = a
    if min_a > a:
        min_a = a
    if max_e < ecc:
        max_e = ecc
    if min_e > ecc:
        min_e = ecc

print ("Processed %d asteroids, a=%f..%f [AU], eccentricity = %f..%f" % (cnt, min_a, max_a, min_e, max_e))

Processed 1024421 asteroids, a=0.555412..3549.257149 [AU], eccentricity = 0.000000..0.998747
In [4]:
# This is the major processing part (that sadly takes a lot of time to complete). It's possible to limit
# the number of asteroids processed by setting max to some value e.g. 20000. However, this will process
# the first 20000 lines from the input file and they're in random (or historical detection) order.
#
# The result is stored in dv_tbl array.
max = 2000000
t0 = time.time()
i = 0

dv_tbl = []
print("Starting calculation, time %s" % time.time())
for e in elements:
    dv_burn1, dv_burn2, dv, tof = calc_min_delta_v(1.0*u.au, e)
    dv_tbl.append([e[0], dv, tof, dv_burn1, dv_burn2])
    i = i + 1
    if (i % 10000 == 0):
        print("Calculated %d orbits, time %s" % (i, time.time()))
    if i >= max:
        break
t1 = time.time()

print("Calculated for %d orbits, calculations took %s s." % (i, (t1-t0)))

Starting calculation, time 1605138454.17286
Calculated 10000 orbits, time 1605138461.25387
Calculated 20000 orbits, time 1605138468.2517946
Calculated 30000 orbits, time 1605138475.2271688
Calculated 40000 orbits, time 1605138482.2088675
Calculated 50000 orbits, time 1605138489.205532
Calculated 60000 orbits, time 1605138496.2299635
Calculated 70000 orbits, time 1605138504.975626
Calculated 80000 orbits, time 1605138512.0324328
Calculated 90000 orbits, time 1605138519.0578794
Calculated 100000 orbits, time 1605138526.0392518
Calculated 110000 orbits, time 1605138533.0169706
Calculated 120000 orbits, time 1605138540.2507865
Calculated 130000 orbits, time 1605138547.542095
Calculated 140000 orbits, time 1605138554.78781
Calculated 150000 orbits, time 1605138561.9948137
Calculated 160000 orbits, time 1605138569.1886222
Calculated 170000 orbits, time 1605138576.3872247
Calculated 180000 orbits, time 1605138583.586671
Calculated 190000 orbits, time 1605138590.7797406
Calculated 200000 orbits, time 1605138597.9799204
Calculated 210000 orbits, time 1605138605.210951
Calculated 220000 orbits, time 1605138612.4713416
Calculated 230000 orbits, time 1605138619.7107937
Calculated 240000 orbits, time 1605138626.92364
Calculated 250000 orbits, time 1605138634.148263
Calculated 260000 orbits, time 1605138641.3619924
Calculated 270000 orbits, time 1605138648.601431
Calculated 280000 orbits, time 1605138655.8200922
Calculated 290000 orbits, time 1605138663.0333924
Calculated 300000 orbits, time 1605138670.256894
Calculated 310000 orbits, time 1605138677.4822667
Calculated 320000 orbits, time 1605138686.8824773
Calculated 330000 orbits, time 1605138694.0940194
Calculated 340000 orbits, time 1605138701.2802842
Calculated 350000 orbits, time 1605138708.5334451
Calculated 360000 orbits, time 1605138715.72933
Calculated 370000 orbits, time 1605138722.9357142
Calculated 380000 orbits, time 1605138730.147699
Calculated 390000 orbits, time 1605138737.3335078
Calculated 400000 orbits, time 1605138744.5789945
Calculated 410000 orbits, time 1605138751.8021677
Calculated 420000 orbits, time 1605138759.0124354
Calculated 430000 orbits, time 1605138766.2160065
Calculated 440000 orbits, time 1605138773.4180093
Calculated 450000 orbits, time 1605138780.6138268
Calculated 460000 orbits, time 1605138787.8240244
Calculated 470000 orbits, time 1605138795.0280218
Calculated 480000 orbits, time 1605138802.223684
Calculated 490000 orbits, time 1605138809.4315357
Calculated 500000 orbits, time 1605138816.658817
Calculated 510000 orbits, time 1605138823.86944
Calculated 520000 orbits, time 1605138831.0564067
Calculated 530000 orbits, time 1605138838.2559657
Calculated 540000 orbits, time 1605138845.4547093
Calculated 550000 orbits, time 1605138852.6610904
Calculated 560000 orbits, time 1605138859.8512597
Calculated 570000 orbits, time 1605138867.0502129
Calculated 580000 orbits, time 1605138874.252377
Calculated 590000 orbits, time 1605138881.456473
Calculated 600000 orbits, time 1605138888.6697218
Calculated 610000 orbits, time 1605138895.9118278
Calculated 620000 orbits, time 1605138903.1254318
Calculated 630000 orbits, time 1605138910.3187733
Calculated 640000 orbits, time 1605138920.301091
Calculated 650000 orbits, time 1605138927.5219295
Calculated 660000 orbits, time 1605138934.7250617
Calculated 670000 orbits, time 1605138941.9334705
Calculated 680000 orbits, time 1605138949.1382103
Calculated 690000 orbits, time 1605138956.384035
Calculated 700000 orbits, time 1605138963.6368895
Calculated 710000 orbits, time 1605138970.8743956
Calculated 720000 orbits, time 1605138978.1090624
Calculated 730000 orbits, time 1605138985.3668282
Calculated 740000 orbits, time 1605138992.6015205
Calculated 750000 orbits, time 1605138999.84376
Calculated 760000 orbits, time 1605139007.109898
Calculated 770000 orbits, time 1605139014.3579655
Calculated 780000 orbits, time 1605139021.603315
Calculated 790000 orbits, time 1605139028.8557277
Calculated 800000 orbits, time 1605139036.1336622
Calculated 810000 orbits, time 1605139043.3490834
Calculated 820000 orbits, time 1605139050.5663166
Calculated 830000 orbits, time 1605139057.807478
Calculated 840000 orbits, time 1605139065.0367649
Calculated 850000 orbits, time 1605139072.2335854
Calculated 860000 orbits, time 1605139079.439212
Calculated 870000 orbits, time 1605139086.628264
Calculated 880000 orbits, time 1605139093.8157477
Calculated 890000 orbits, time 1605139101.0100577
Calculated 900000 orbits, time 1605139108.1839676
Calculated 910000 orbits, time 1605139115.3756688
Calculated 920000 orbits, time 1605139122.5915153
Calculated 930000 orbits, time 1605139129.8155963
Calculated 940000 orbits, time 1605139137.035852
Calculated 950000 orbits, time 1605139144.2716377
Calculated 960000 orbits, time 1605139151.5141861
Calculated 970000 orbits, time 1605139158.7397847
Calculated 980000 orbits, time 1605139165.934174
Calculated 990000 orbits, time 1605139173.1366704
Calculated 1000000 orbits, time 1605139180.3339329
Calculated 1010000 orbits, time 1605139187.5351775
Calculated 1020000 orbits, time 1605139194.7177832
Calculated for 1028453 orbits, calculations took 746.6426937580109 s.
In [5]:
# Now sort the table by the second (index 1) value. This is the total delta-v
t2 = time.time()
dv_tbl2 = sorted(dv_tbl, key=lambda e: e[1].value)
t3 = time.time()

print("dv_tbl (%d elements) is now sorted, sorting took %s s." % (len(dv_tbl), t3-t2))

dv_tbl (1028453 elements) is now sorted, sorting took 1.0772876739501953 s.
In [6]:
# Now print the top 10 asteroids (with smallest delta-v)
top_x = 10
print(len(dv_tbl2))
print(" TOP %d ASTEROIDS THAT ARE EASIEST TO REACH" % top_x)
for i in range(top_x):
    x = dv_tbl2[i]
    print("%d. %s:  fly-by delta-v=%s randezvous delta-v=%s" % (i+1, x[0], x[3], x[1] ))

1028453
 TOP 10 ASTEROIDS THAT ARE EASIEST TO REACH
1. 2010 FN:  fly-by delta-v=100.73845547422833 m / s randezvous delta-v=178.38592193633747 m / s
2. 2018 UC:  fly-by delta-v=132.27082396382207 m / s randezvous delta-v=189.4751242873342 m / s
3. 2016 JA:  fly-by delta-v=111.15564339177044 m / s randezvous delta-v=218.82623504831835 m / s
4. 2020 SO:  fly-by delta-v=123.50356470013689 m / s randezvous delta-v=224.53306310902377 m / s
5. 2019 YB4:  fly-by delta-v=280.1225044691703 m / s randezvous delta-v=308.8324268322133 m / s
6. 2019 AU6:  fly-by delta-v=218.92547156701437 m / s randezvous delta-v=342.59523498058877 m / s
7. 2000 SG344:  fly-by delta-v=180.04330838593847 m / s randezvous delta-v=351.4091885944977 m / s
8. 2010 JK1:  fly-by delta-v=232.52398230729133 m / s randezvous delta-v=446.5553471761266 m / s
9. 2016 WQ3:  fly-by delta-v=348.1024745741976 m / s randezvous delta-v=466.27559915444874 m / s
10. 2013 BS45:  fly-by delta-v=406.0127279657208 m / s randezvous delta-v=468.60941121661375 m / s
In [7]:
# And write all the data to output file.
f = open("../data/mpc/asteroids.txt", "w+")
f.write("#name, total delta-v[m/s], time of flight[days], Hohmann burn1 [m/s], Hohmann burn2 [m/s]\n")
for x in dv_tbl2:
    f.write("%s, %4.2f, %4.1f, %4.2f, %4.2f\n" % (x[0], x[1].value, x[2].value, x[3].value, x[4].value))
f.close()