Kepler orbital elements with plotly

This was an experiment to see how hard it was to create an interactive figure with plotly.

The idea is to adjust an orbit's Keplerian elements using sliders to get a quick feel for what they mean. Static pictures can get the point across, but playing with the elements wil lbe the quickest way to remember what they mean when learning for the first time.

To run the widget, click the "Run" button on the cell below.

To make the slider widget, I followed nearly the same format as the ipywidget example: https://ipywidgets.readthedocs.io/en/latest/examples/Using%20Interact.html and call my kep2cart on the output of sliders to convert them to x, y, z that I plot for one orbit.



In [10]:

    
def rot(angle, axis, in_degrees=False):
    """
    Find a 3x3 euler rotation matrix given an angle and axis.

    Rotation matrix used for rotating a vector about a single axis.

    Args:
        angle (float): angle in degrees to rotate
        axis (int): 1, 2 or 3
        in_degrees (bool): specify the angle in degrees. if false, using
            radians for `angle`
    """
    R = np.eye(3)
    if in_degrees:
        angle = np.deg2rad(angle)
    cang = cos(angle)
    sang = sin(angle)
    if axis == 1:
        R[1, 1] = cang
        R[2, 2] = cang
        R[1, 2] = sang
        R[2, 1] = -sang
    elif axis == 2:
        R[0, 0] = cang
        R[2, 2] = cang
        R[0, 2] = -sang
        R[2, 0] = sang
    elif axis == 3:
        R[0, 0] = cang
        R[1, 1] = cang
        R[1, 0] = -sang
        R[0, 1] = sang
    else:
        raise ValueError("axis must be 1, 2 or 2")
    return R


def R1(angle, in_degrees=False):
    return rot(angle, 1, in_degrees=in_degrees)

def R3(angle, in_degrees=False):
    return rot(angle, 3, in_degrees=in_degrees)


def kep2cart(mu, kep_oe_vector):
    # Steps outline:
    # 1. Using sma, ecc, tru, solve for r, v in the PQW frame
    # 2. Check singularities, make transform matrix T_PQW_ijk
    # 3. transform r and v fromPQW to inertial

    # 1. Using sma, ecc, tru, solve for r, v in the PQW frame
    sma, ecc, inc, raan, argp, tru = kep_oe_vector

    p = sma * (1 - ecc**2)  # Semiparameter

    rPQW = np.array([p * cos(tru), p * sin(tru), 0]) / (1 + ecc * cos(tru))
    vPQW = np.array([-sin(tru), ecc + cos(tru), 0]) * sqrt(mu / p)

    # 2. Check singularities, make transform matrix T_PQW_ijk
    # with raan, inc, argp using 3-1-3 Euler sequence
    if ecc < 1e-12 and inc < 1e-12:
        # raan = 0; argp = 0; tru = l;
        T_pqw_ijk = np.eye(3)
    elif ecc < 1e-12:
        # argp = 0; tru = u;
        T_pqw_ijk = R3(-raan) @ R1(-inc)
    elif inc < 1e-12:
        # raan = 0,  argp = omegaBar;
        T_pqw_ijk = R3(-argp)
    else:
        T_pqw_ijk = R3(-raan) @ R1(-inc) @ R3(-argp)

    # 3. transform r and v fromPQW to inertial
    r_vec = T_pqw_ijk @ rPQW
    v_vec = T_pqw_ijk @ vPQW

    return r_vec, v_vec

max_lim = 10

scene=dict(camera=dict(eye=dict(x=1.15, y=1.15, z=0.8)), #the default values are 1.25, 1.25, 1.25
   xaxis={'range': [-max_lim, max_lim]},
   yaxis={'range': [-max_lim, max_lim]},
   zaxis={'range': [-max_lim, max_lim]},
   aspectmode='cube', #this string can be 'data', 'cube', 'auto', 'manual'
   #a custom aspectratio is defined as follows:
   aspectratio=dict(x=1, y=1, z=0.95)
)

# eye gives the position of the camera eye;
# aspectmode='cube', the scene’s axes are drawn as a cube, regardless of the axes’ ranges
# aspectmode='data' preserves the proportion of axes ranges
# 'manual' when you set the aspectratio,
# 'auto' the scene’s axes are drawn with 'data', except when one axis is more than four times
# the size of the two others; in that case the 'cube' is used.

fig = go.FigureWidget(**{'layout': dict(
    scene=scene,
    width=700,
    height=700,
    )}
)

scatt = fig.add_scatter3d()
scatt.marker.size = 2
scatt.line.color = 'blue'

mu = 1

sma_tup = (1, 5, .3)  # min, max, step
ecc_tup = (0, .8, .1)
inc_tup = (0, np.pi/2, .1)  
raan_tup = (0, np.pi, .1)  
argp_tup = (0, np.pi, .1)  

tru_array = np.linspace(0, 2*np.pi, 50)
cur_r_vecs = np.empty((3, len(tru_array)))

@interact(sma=sma_tup, ecc=ecc_tup, inc=inc_tup, raan=raan_tup, argp=argp_tup)
def update(sma=1., ecc=0., inc=0., raan=0.0, argp=0.0,):
    with fig.batch_update():
        for idx, tru in enumerate(tru_array):
            r_vec, v_vec = kep2cart(mu, np.array([sma, ecc, inc, raan, argp, tru]))
            cur_r_vecs[:, idx] = r_vec
            
        xs, ys, zs = cur_r_vecs
        scatt.x=xs
        scatt.y=ys
        scatt.z=zs

fig

Kepler orbital elements with plotly

The plotly widget