Python EB

All Planets

import bpy
import bmesh
import math
import numpy

bpyscene = bpy.context.scene

#___________________Sun______________________

# empty mesh and the planet
mesh = bpy.data.meshes.new('Sun_Sphere')
sun_sphere = bpy.data.objects.new("Sun_Sphere", mesh)

# adding the planet
bpy.context.collection.objects.link(sun_sphere)
bpy.context.view_layer.objects.active = sun_sphere
bpy.context.active_object.select_set(state=True)

# constructing bmesh sphere and assigning it to the blender mesh
bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.1)
bm.to_mesh(mesh)
bm.free()

# changing mesh properties - ?
bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

# making an empty object to represent a reference frame of the planet
def MakeEmpty(Name):
    e = None
    for i in bpy.data.objects:
        if i.name == Name: e = i
    if e == None:
        bpy.ops.object.add(type='EMPTY')
        e = bpy.context.active_object
        e.name = Name
    return(e)

MakeEmpty('Sun')

# making the reference frame parent of the planet to stick them together
objects = bpy.data.objects
a1 = objects['Sun_Sphere']
b1 = objects['Sun']
a1.parent = b1

periloc = (0,0,0)

# transforming the location and rotation of the reference frame; placing the planet at perihelion by default
b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-7.2*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-346.6*math.pi/180,orient_axis='Z')


#___________________Mercury_____________________

mesh = bpy.data.meshes.new('Mercury_Sphere')
mercury_sphere = bpy.data.objects.new("Mercury_Sphere", mesh)

bpy.context.collection.objects.link(mercury_sphere)
bpy.context.view_layer.objects.active = mercury_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Mercury')

objects = bpy.data.objects
a1 = objects['Mercury_Sphere']
b1 = objects['Mercury']
a1.parent = b1

# calculating the location of the perihelion to place the planet there
a = 0.3870978
ecc = 0.2056386
inc = 7.00393*math.pi/180
asc = 48.3082*math.pi/180
peri = 29.1753*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-7.1*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-318.8*math.pi/180,orient_axis='Z')

#____________________Venus_____________________

mesh = bpy.data.meshes.new('Venus_Sphere')
venus_sphere = bpy.data.objects.new("Venus_Sphere", mesh)

bpy.context.collection.objects.link(venus_sphere)
bpy.context.view_layer.objects.active = venus_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Venus')

objects = bpy.data.objects
a1 = objects['Venus_Sphere']
b1 = objects['Venus']
a1.parent = b1

a = 0.7233304
ecc = 0.006802
inc = 3.39449*math.pi/180
asc = 76.62845*math.pi/180
peri = 54.76321*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-178.8*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-210.18*math.pi/180,orient_axis='Z')

#____________________Earth_____________________

mesh = bpy.data.meshes.new('Earth_Sphere')
earth_sphere = bpy.data.objects.new("Earth_Sphere", mesh)

bpy.context.collection.objects.link(earth_sphere)
bpy.context.view_layer.objects.active = earth_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Earth')

objects = bpy.data.objects
a1 = objects['Earth_Sphere']
b1 = objects['Earth']
a1.parent = b1

a = 1
ecc = 0.0167
inc = 0*math.pi/180
asc = 174.4519*math.pi/180
peri = 288.1*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-23.4*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-90*math.pi/180,orient_axis='Z')

#____________________Mars_____________________

mesh = bpy.data.meshes.new('Mars_Sphere')
mars_sphere = bpy.data.objects.new("Mars_Sphere", mesh)

bpy.context.collection.objects.link(mars_sphere)
bpy.context.view_layer.objects.active = mars_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Mars')

objects = bpy.data.objects
a1 = objects['Mars_Sphere']
b1 = objects['Mars']
a1.parent = b1

a = 1.5237
ecc = 0.093417
inc = 1.84833*math.pi/180
asc = 49.508*math.pi/180
peri = 286.60205*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-26.7*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-353.6*math.pi/180,orient_axis='Z')

#____________________Ceres_____________________

mesh = bpy.data.meshes.new('Ceres_Sphere')
ceres_sphere = bpy.data.objects.new("Ceres_Sphere", mesh)

bpy.context.collection.objects.link(ceres_sphere)
bpy.context.view_layer.objects.active = ceres_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Ceres')

objects = bpy.data.objects
a1 = objects['Ceres_Sphere']
b1 = objects['Ceres']
a1.parent = b1

a = 2.7671
ecc = 0.07610
inc = 10.60069*math.pi/180
asc = 80.658515*math.pi/180
peri = 71.44921*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-8.4*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-12.1*math.pi/180,orient_axis='Z')

#____________________Jupiter_____________________

mesh = bpy.data.meshes.new('Jupiter_Sphere')
jupiter_sphere = bpy.data.objects.new("Jupiter_Sphere", mesh)

bpy.context.collection.objects.link(jupiter_sphere)
bpy.context.view_layer.objects.active = jupiter_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Jupiter')

objects = bpy.data.objects
a1 = objects['Jupiter_Sphere']
b1 = objects['Jupiter']
a1.parent = b1

a = 5.20336
ecc = 0.04839
inc = 1.3053*math.pi/180
asc = 100.55615*math.pi/180
peri = 273.62*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-2.2*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-248.9*math.pi/180,orient_axis='Z')

#____________________Saturn_____________________

mesh = bpy.data.meshes.new('Saturn_Sphere')
saturn_sphere = bpy.data.objects.new("Saturn_Sphere", mesh)

bpy.context.collection.objects.link(saturn_sphere)
bpy.context.view_layer.objects.active = saturn_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Saturn')

objects = bpy.data.objects
a1 = objects['Saturn_Sphere']
b1 = objects['Saturn']
a1.parent = b1

a = 9.53707
ecc = 0.05415
inc = 2.48446*math.pi/180
asc = 113.715*math.pi/180
peri = 340.0014*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-28.1*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-80.2*math.pi/180,orient_axis='Z')

#____________________Uranus_____________________

mesh = bpy.data.meshes.new('Uranus_Sphere')
uranus_sphere = bpy.data.objects.new("Uranus_Sphere", mesh)

bpy.context.collection.objects.link(uranus_sphere)
bpy.context.view_layer.objects.active = uranus_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Uranus')

objects = bpy.data.objects
a1 = objects['Uranus_Sphere']
b1 = objects['Uranus']
a1.parent = b1

a = 19.1323
ecc = 0.049295
inc = 0.76986*math.pi/180
asc = 74.2298*math.pi/180
peri = 99.3024*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-82.3*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-258.3*math.pi/180,orient_axis='Z')

#____________________Neptune_____________________

mesh = bpy.data.meshes.new('Neptune_Sphere')
neptune_sphere = bpy.data.objects.new("Neptune_Sphere", mesh)

bpy.context.collection.objects.link(neptune_sphere)
bpy.context.view_layer.objects.active = neptune_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Neptune')

objects = bpy.data.objects
a1 = objects['Neptune_Sphere']
b1 = objects['Neptune']
a1.parent = b1

a = 30.0485
ecc = 0.0067291
inc = 1.77759*math.pi/180
asc = 131.9249*math.pi/180
peri = 269.5025*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-28.4*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-318.6*math.pi/180,orient_axis='Z')

#____________________Pluto_____________________

mesh = bpy.data.meshes.new('Pluto_Sphere')
pluto_sphere = bpy.data.objects.new("Pluto_Sphere", mesh)

bpy.context.collection.objects.link(pluto_sphere)
bpy.context.view_layer.objects.active = pluto_sphere
bpy.context.active_object.select_set(state=True)

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm, u_segments=64, v_segments=64, diameter=0.05)
bm.to_mesh(mesh)
bm.free()

bpy.ops.object.modifier_add(type='SUBSURF')
bpy.ops.object.shade_smooth()

MakeEmpty('Pluto')

objects = bpy.data.objects
a1 = objects['Pluto_Sphere']
b1 = objects['Pluto']
a1.parent = b1

a = 39.4816
ecc = 0.2488
inc = 17.1335*math.pi/180
asc = 110.3108*math.pi/180
peri = 114.5568*math.pi/180
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

periloc = numpy.dot(Eul,(a*(1 - ecc), 0, 0))

b1.location = periloc
bpy.context.active_object.select_set(state=True)
bpy.ops.transform.rotate(value=-112.8*math.pi/180,orient_axis='Y')
bpy.ops.transform.rotate(value=-138.1*math.pi/180,orient_axis='Z')

Coordinate Systems

import bpy
import math
from mathutils import Matrix

# this script generates three empties, one of which is the galactic coordinate system, 
# the other is equatorial, and the third is HGI (solar);
# thus we might draw in any of these and then just parent to the proper empty, 
# and it would show up properly in our default ecliptic system

#______________Galactic coordinate system______________

# inclination of the galactic plane
inc = 60.2*math.pi/180

# longitude of the ascending node
asc = 270.1*math.pi/180

# argument of X axis - towards the Sun's motion in the Galaxy 
peri = 353.5*math.pi/180

# making an empty object to represent the galactic reference frame
def MakeEmpty(Name):
    e = None
    for i in bpy.data.objects:
        if i.name == Name: e = i
    if e == None:
        bpy.ops.object.add(type='EMPTY')
        e = bpy.context.active_object
        e.name = Name
    return(e)

MakeEmpty('Galactic_system')

# activating the exmpty and rotating it
bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(asc, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(inc, 3, 'X')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(peri, 3, 'Z')).to_euler()

#______________Ecliptic coordinate system______________

inc = 23.4*math.pi/180
asc = 180*math.pi/180
peri = 180*math.pi/180

MakeEmpty('Ecliptic_system')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(asc, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(inc, 3, 'X')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(peri, 3, 'Z')).to_euler()

#______________HGI coordinate system______________

inc = 7.25*math.pi/180
asc = 75.767*math.pi/180
peri = 0*math.pi/180

MakeEmpty('HGI_system')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(asc, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(inc, 3, 'X')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(peri, 3, 'Z')).to_euler()

Direction Arrows

import bpy
import math
import bmesh
import mathutils
from mathutils import Matrix
from mathutils import Vector

# generating some important directions, by default they're linked to the appropriate object's Z axis

# a function that gives you the arrow
def gimmearrow(Arrowname):
    for ob in bpy.context.selected_objects:
        bpy.context.active_object.select_set(state=False)
    obj = bpy.context.window.scene.objects[0]
    bpy.context.view_layer.objects.active = None
    # start with a cone and a cylinder
    bpy.ops.mesh.primitive_cone_add(vertices=32, radius1=0.03, radius2=0.0, depth=0.06, location=(0.0, 0.0, 1.0), rotation=(0.0, 0.0, 0.0))
    for obj in bpy.context.selected_objects:
        obj.name = "con"
    cyl = bpy.ops.mesh.primitive_cylinder_add(vertices=32, radius=0.012, depth=1.0,location=(0.0, 0.0, 0.5), rotation=(0.0, 0.0, 0.0))
    for obj in bpy.context.selected_objects:
        obj.name = "cyl"
    # selecting the bastards and merging them    
    bpy.data.objects['con'].select_set(True)
    bpy.data.objects['cyl'].select_set(True)
    bpy.ops.object.join()
    # renaming the joined object to the input name
    for obj in bpy.context.selected_objects:
        obj.name = Arrowname
        obj.data.name = "Arrow"

# the infamous makeempty function
def MakeEmpty(Name):
    e = None
    for i in bpy.data.objects:
        if i.name == Name: e = i
    if e == None:
        bpy.ops.object.add(type='EMPTY')
        e = bpy.context.active_object
        e.name = Name
    return(e)

# and finally some action

#_____________Solar apex____________
# H. Karttunen et al. Fundamental Astronomy. 2007. P. 350.

gimmearrow("Direction")

lon = 270*math.pi/180
lat = 53.4*math.pi/180

MakeEmpty('Solar_apex')

# activating the empty and rotating it
bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

#parenting
objects = bpy.data.objects
a1 = objects['Direction']
b1 = objects['Solar_apex']
a1.parent = b1

#_____________Interstellar field____________
# E. J. Zirnstein et al. Local Interstellar Magnetic Field Determined From the Interstellar Boundary Explorer Ribbon. APJ Letters, 818, No 1. 2016.
# alternatively - IceCube (https://arxiv.org/abs/1812.05682) gives lon 223.2, lat 28.5
# and light polarization maps - https://arxiv.org/abs/1806.02806 - give lon 247 lat 82.2
# and Cassini (see DOI: 10.1088/0004-637X/778/1/40) gives lon 190 lat 15 
  
gimmearrow("Direction2")

lon = 227.28*math.pi/180
lat = 34.62*math.pi/180

MakeEmpty('Interstellar_mag_field')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction2']
b1 = objects['Interstellar_mag_field']
a1.parent = b1

#____________Cassini/INCA belt center____________
# S. M. Krimigis et al. Imaging the Interaction of the Heliosphere with the Interstellar Medium from Saturn with Cassini. Science, 326, 2009, P. 971-973

gimmearrow("Direction3")

lon = 15*math.pi/180
lat = -10*math.pi/180

MakeEmpty('Cassini_belt')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction3']
b1 = objects['Cassini_belt']
a1.parent = b1

#___________Heliospheric nose_______________
# https://arxiv.org/pdf/1111.3075.pdf

gimmearrow("Direction4")

lon = 255*math.pi/180
lat = 5*math.pi/180

MakeEmpty('Heliospheric_nose')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction4']
b1 = objects['Heliospheric_nose']
a1.parent = b1

#___________Interstellar wind_______________
# https://arxiv.org/pdf/1501.02725.pdf

gimmearrow("Direction5")

lon = 75*math.pi/180
lat = -5*math.pi/180

MakeEmpty('Interstellar_wind')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction5']
b1 = objects['Interstellar_wind']
a1.parent = b1

#____________IBEX ribbon center____________
# H. O. Funsten et al. Structures and Spectral Variations of the Outer Heliosphere in IBEX Energetic Neutral Atom Maps

gimmearrow("Direction6")

lon = 221*math.pi/180
lat = 39*math.pi/180

MakeEmpty('IBEX_ribbon')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction6']
b1 = objects['IBEX_ribbon']
a1.parent = b1

#____________Fine-structure hot pole____________
# J. K. Webb et al. Indications of a Spatial Variation of the Fine Structure Constant
# though https://arxiv.org/pdf/1412.3418.pdf gives lon 248 lat -57.7

gimmearrow("Direction7")

lon = 265.2*math.pi/180
lat = -34.6*math.pi/180

MakeEmpty('Alpha_pole')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction7']
b1 = objects['Alpha_pole']
a1.parent = b1

#____________CMB hot pole____________
# http://www.astro.ucla.edu/~wright/CMB-DT.html

gimmearrow("Direction8")

lon = 171.6*math.pi/180
lat = -11.1*math.pi/180

MakeEmpty('CMB_dipole')

bpy.context.active_object.select_set(state=True)
obj = bpy.context.active_object

obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(lon, 3, 'Z')).to_euler()
obj.rotation_euler = (obj.rotation_euler.to_matrix() @ Matrix.Rotation(90*math.pi/180-lat, 3, 'Y')).to_euler()

objects = bpy.data.objects
a1 = objects['Direction8']
b1 = objects['CMB_dipole']
a1.parent = b1

Orbits

import bpy
import math
import numpy
from mathutils import Vector  

#________________ Mercury__________________

# semimajor axis
a = 0.3870978

# eccentricity
ecc = 0.2056386

# inclination
inc = 7.00393*math.pi/180

# longitude of the ascending node
asc = 48.3082*math.pi/180

# argument of the perihelion
peri = 29.1753*math.pi/180

# Euler matrix to rotate the frame so the orbit would lie in [the new] XoY plane
Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

# number of points to generate the curve
steps = 100

# generating the "flat" ellipse, not yet with the proper orientation
plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

# generation of the points - applying Euler matrix to the ellipse equation
listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

# ^ so tuple transforms [] type of array into () (to be used by Vector), 
# ndarray.tolist eliminates some of the internal nesting, 
# and numpy.dot is just a matrix multiplication

# weight (for the curve points) 
w = 1 
 
# defines the procedure to draw a curve 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True
    
# creates the curve
MakePolyLine("Mercury", "Orbit", listOfVectors)

#__________________Venus____________________

a = 0.7233304
ecc = 0.006802
inc = 3.39449*math.pi/180
asc = 76.62845*math.pi/180
peri = 54.76321*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Venus", "Orbit", listOfVectors)

#__________________Earth_____________________

a = 1
ecc = 0.0167
inc = 0*math.pi/180
asc = 174.4519*math.pi/180
peri = 288.1*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Earth", "Orbit", listOfVectors)

#___________________Mars____________________

a = 1.5237
ecc = 0.093417
inc = 1.84833*math.pi/180
asc = 49.508*math.pi/180
peri = 286.60205*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Mars", "Orbit", listOfVectors)

#___________________Ceres____________________

a = 2.7671
ecc = 0.07610
inc = 10.60069*math.pi/180
asc = 80.658515*math.pi/180
peri = 71.44921*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Ceres", "Orbit", listOfVectors)


#___________________Jupiter____________________

a = 5.20336
ecc = 0.04839
inc = 1.3053*math.pi/180
asc = 100.55615*math.pi/180
peri = 273.62*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Jupiter", "Orbit", listOfVectors)

#___________________Saturn____________________

a = 9.53707
ecc = 0.05415
inc = 2.48446*math.pi/180
asc = 113.715*math.pi/180
peri = 340.0014*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Saturn", "Orbit", listOfVectors)

#___________________Uranus____________________

a = 19.1323
ecc = 0.049295
inc = 0.76986*math.pi/180
asc = 74.2298*math.pi/180
peri = 99.3024*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Uranus", "Orbit", listOfVectors)

#___________________Neptune____________________

a = 30.0485
ecc = 0.0067291
inc = 1.77759*math.pi/180
asc = 131.9249*math.pi/180
peri = 269.5025*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Neptune", "Orbit", listOfVectors)

#___________________Pluto____________________

a = 39.4816
ecc = 0.2488
inc = 17.1335*math.pi/180
asc = 110.3108*math.pi/180
peri = 114.5568*math.pi/180

Eul = numpy.array([[math.cos(asc)*math.cos(peri) - math.cos(inc)*math.sin(asc)*math.sin(peri),-math.cos(inc)*math.cos(peri)*math.sin(asc) - math.cos(asc)
*math.sin(peri), math.sin(asc)*math.sin(inc)],
[math.cos(peri)*math.sin(asc) + math.cos(asc)*math.cos(inc)*math.sin(peri), math.cos(asc)*math.cos(inc)*math.cos(peri) - math.sin(asc)*math.sin(peri), -
math.cos(asc)*math.sin(inc)],
[math.sin(inc)*math.sin(peri), math.cos(peri)*math.sin(inc), math.cos(inc)]])

steps = 100

plist = [None] * steps
i = 0
while i < steps:
    plist[i] = [a*(math.cos(2*math.pi*i/(steps-1)) - ecc), a*math.sqrt(1 - ecc**2)*math.sin(2*math.pi*i/(steps-1)), 0]
    i += 1

listOfVectors = [None] * steps
i = 0
while i < steps:
    listOfVectors[i] = tuple(numpy.ndarray.tolist(numpy.dot(Eul,plist[i])))
    i += 1

w = 1 
 
def MakePolyLine(objname, curvename, cList):  
    curvedata = bpy.data.curves.new(name=curvename, type='CURVE')  
    curvedata.dimensions = '3D'  
  
    objectdata = bpy.data.objects.new(objname, curvedata)  
    objectdata.location = (0,0,0) #object origin  
    bpy.context.collection.objects.link(objectdata)  
  
    polyline = curvedata.splines.new('NURBS')  
    polyline.points.add(len(cList)-1)  
    for num in range(len(cList)):  
        polyline.points[num].co = (cList[num])+(w,)  
  
    polyline.order_u = len(polyline.points)-1
    polyline.use_endpoint_u = True

MakePolyLine("Pluto", "Orbit", listOfVectors)