gistfile1.txt

Traceback (most recent call last):
  File "C:/Users/taiko/PycharmProjects/untitled/solar_system_toy.py", line 125, in <module>
    orbit_list[idx] = orbit_list[idx].propagate(1*u.day,method=cowell,ad=ad)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\poliastro\twobody\orbit.py", line 999, in propagate
    cartesian = propagate(self, time_of_flight, method=method, rtol=rtol, **kwargs)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\poliastro\twobody\propagation.py", line 468, in propagate
    orbit.attractor.k, orbit.r, orbit.v, time_of_flight.to(u.s), rtol=rtol, **kwargs
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\poliastro\twobody\propagation.py", line 105, in cowell
    dense_output=True,
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 477, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\poliastro\integrators.py", line 375, in __init__
    self.f = self.fun(self.t, self.y)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\integrate\_ivp\base.py", line 139, in fun
    return self.fun_single(t, y)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\integrate\_ivp\base.py", line 21, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)
  File "C:\Users\taiko\AppData\Local\Programs\Python\Python37-32\lib\site-packages\poliastro\core\propagation.py", line 37, in func_twobody
    ax, ay, az = ad(t0, u_, k, **ad_kwargs)
TypeError: 'Quantity' object is not callable

solar_system_toy.py

from astropy.time import Time
from poliastro.twobody import Orbit
from poliastro.bodies import  Sun
from astropy import units as u
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from skyfield.api import  load
from scipy.spatial import distance
from poliastro.twobody.propagation import propagate, cowell

def solar_pertubation(planet_orbit_list, av):
    m_list = [3.301e23*u.kg, 4.867e24*u.kg, 5.973e24*u.kg, 6.417e23*u.kg, 1.899e27*u.kg, 5.685e26*u.kg, 8.682e25*u.kg,
              1.024e26*u.kg]  # -> mass of planet 1 , mass of planet 2 .... units

    total_accel = 0

    for i in range(len(planet_orbit_list)):
        if (planet_orbit_list[i].r.value == av.r.value).all():
            return None
        else:
            delta_r = planet_orbit_list[i].r - av.r

            total_accel += (G * m_list[i] * delta_r / (distance.euclidean(planet_orbit_list[i].r, av.r)) ** 3)#.value

    return total_accel

G = 6.67408e-11 #*u.m**3/(u.kg*u.s**2)#.to(u.km ** 3 / u.s ** 2)#untis

solar_system_r = []
solar_system_v = []

sun_r = []
sun_v = []

ts = load.timescale()
t = ts.tai_jd(2458689.5)  #ts.now() #
epoch = Time("2019-07-25 12:05:50", scale="tdb")

planets = load('de421.bsp')
#sun= planets['sun']
mercury = planets['mercury']
venus = planets['venus']
earth = planets['earth']
mars = planets['mars']
jupiter = planets['JUPITER BARYCENTER']
saturn = planets['saturn BARYCENTER']
uranus = planets['uranus BARYCENTER']
neptune = planets['neptune BARYCENTER']
pluto = planets['pluto BARYCENTER']

#sun_r.append(sun.at(t).position.km)
solar_system_r.append(mercury.at(t).position.km)
solar_system_r.append(venus.at(t).position.km)
solar_system_r.append(earth.at(t).position.km)
solar_system_r.append(mars.at(t).position.km)
solar_system_r.append(jupiter.at(t).position.km)
solar_system_r.append(saturn.at(t).position.km)
solar_system_r.append(uranus.at(t).position.km)
solar_system_r.append(neptune.at(t).position.km)


#sun_v.append(sun.at(t).velocity.km_per_s)
solar_system_v.append(mercury.at(t).velocity.km_per_s)
solar_system_v.append(venus.at(t).velocity.km_per_s)
solar_system_v.append(earth.at(t).velocity.km_per_s)
solar_system_v.append(mars.at(t).velocity.km_per_s)
solar_system_v.append(jupiter.at(t).velocity.km_per_s)
solar_system_v.append(saturn.at(t).velocity.km_per_s)
solar_system_v.append(uranus.at(t).velocity.km_per_s)
solar_system_v.append(neptune.at(t).velocity.km_per_s)


planet_names = ['Mercury', 'Venus','Earth','Mars','Jupiter','Saturn','Uranus','Neptune','Pluto'] #'Sun',
stars = ['Sun']

a = np.load('asteroid_data_new.npy', allow_pickle=True).tolist()
df = pd.read_csv('file_name.csv',header=None)
df[['vx','vy','vz','rx','ry','rz']] = df[0].str.split(' ', expand=True)
df = df.drop(columns=[0])


df_v = np.array(df.iloc[:50, :3]) # df.iloc[1:i:3] where i is the number of asteroids you want to plot velocity vectors x,y,z
df_r = np.array(df.iloc[:50, 3:]) # df.iloc[1:i:3] where i is the number of asteroids you want to plot position vectors x,y,z
df_v =df_v.astype(float)
df_r =df_r.astype(float)


epoch = Time("2019-07-25 12:05:50", scale="tdb")
df = pd.DataFrame(a)
epoch_data = df['Epoch']
orbit_list = []
planet_orbit_list = []

for i, j in zip(solar_system_r, solar_system_v):
    i = i * u.km
    j = j * u.km / u.s

    planet_orbit = Orbit.from_vectors(Sun, i, j, epoch=epoch)

    planet_orbit_list.append(planet_orbit)

for i,j in zip(df_r,df_v):

      i = i * u.km
      j = j * u.km/ u.s

      initial = Orbit.from_vectors(Sun, i, j,epoch=epoch)

      orbit_list.append(initial)


n = 0
while n < 1:

    plt.style.use('dark_background')
    axes = plt.gca()
    axes.set_xlim([-1e9, 1e9])  # km axis
    axes.set_ylim([-1e9, 1e9])

    plt.plot(0,0, marker='o', markersize=2,color='yellow')

    for idx in range(len(orbit_list)):
        ad = solar_pertubation(planet_orbit_list,orbit_list[idx])
        orbit_list[idx] = orbit_list[idx].propagate(1*u.day,method=cowell,ad=ad)


        plt.plot(orbit_list[idx].r.to_value().item(0), orbit_list[idx].r.to_value().item(1), marker='o',ms=72./300,markeredgewidth=0, color='green')

    for idx in range(len(planet_orbit_list)):
        planet_orbit_list[idx] = planet_orbit_list[idx].propagate(1 * u.day)
        plt.plot(planet_orbit_list[idx].r.to_value().item(0), planet_orbit_list[idx].r.to_value().item(1), marker='o', markersize=0.5,markeredgewidth=0)
    n += 1
    plt.xlabel('[km]', fontsize=18)
    plt.ylabel('[km]', fontsize=16)
    plt.title(str(orbit_list[1].epoch)[:-12])#remove last 12 chars


    plt.savefig( str(n) + '.png',dpi=300)
    plt.clf()
    print('propagating '+ str(n) + ' days...')