In [1]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
rc('text', usetex=True)
plt.rc('font', family='serif')

from __future__ import print_function


from IPython.html.widgets import interact


    

In [2]:

    

%%bash 
make
./Lindblad_resonances


    

    




gfortran -O2 -c physics_Lindblad_resonances.f08
gfortran -O2 -c ../bulirsch_stoer.f08 
gfortran -O2 -o Lindblad_resonances  physics_Lindblad_resonances.f08 ../bulirsch_stoer.f08 Lindblad_resonances.f08
 energy (initial): 
  -1.6812746142753989     
 starting body: 
           1
 Energy of body           1
 -0.52506248501973873     
 starting body: 
           2
 Energy of body           2
 -0.65871916014014653     
 starting body: 
           3
 Energy of body           3
 -0.49749296911551372     
 
 energy (final): 
  -1.6812746084821135     
 overall conservation:
 energy err:
   3.4457698687426013E-009

In [3]:

    

n_bodies = 3

i=1
t_i, x_i, y_i, z_i, v_x_i, v_y_i, v_z_i = np.loadtxt('body_' + str(i).zfill(4) + '.dat', unpack=True, skiprows=1)
n_points = t_i.size

t = np.empty( (n_bodies, n_points))
t[i-1,:]   = t_i
x = np.empty( (n_bodies, n_points))
x[i-1,:]   = x_i
y = np.empty( (n_bodies, n_points))
y[i-1,:]   = y_i
z = np.empty( (n_bodies, n_points))
z[i-1,:]   = z_i
v_x = np.empty( (n_bodies, n_points))
v_x[i-1,:] = v_x_i
v_y = np.empty( (n_bodies, n_points))
v_y[i-1,:] = v_y_i
v_z = np.empty( (n_bodies, n_points))
v_z[i-1,:] = v_z_i 


for i in range(2, n_bodies+1):
    t_i, x_i, y_i, z_i, v_x_i, v_y_i, v_z_i = np.loadtxt('body_' + str(i).zfill(4) + '.dat', unpack=True, skiprows=1)

    t[i-1,:]   = t_i
    x[i-1,:]   = x_i
    y[i-1,:]   = y_i
    z[i-1,:]   = z_i
    v_x[i-1,:] = v_x_i
    v_y[i-1,:] = v_y_i
    v_z[i-1,:] = v_z_i


    

In [4]:

    

plot_step_size = int(t[0,-1] - t[0, 0]) / 10
def plot_solution(plot_start=0):
    print("plot range: t= ",plot_start,  " - ", plot_start + plot_step_size)
    plt.figure()

    
    indices = np.where( ((t[0,:] > plot_start) & (t[0,:] < plot_start + plot_step_size)))
    plt.subplot(121, aspect="equal")
    plt.plot(x[0][indices], y[0][indices])
    plt.scatter(0,0, marker="+")
    plt.title("Corotational Resonance")
    plt.xlabel(r"$x$ (along bar)")
    plt.ylabel(r"$y$ (perpendicular to bar")
    plt.xlim(-2,2)
    plt.ylim(-2,2)

    
    plt.subplot(122, aspect="equal")
    plt.plot(x[1][indices], y[1][indices])
    plt.scatter(0,0, marker="+")
    plt.title("Outer Lindblad Resonance")
    plt.xlabel(r"$x$ (along bar)")
    plt.ylabel(r"$y$ (perpendicular to bar")
    plt.xlim(-2,2)
    plt.ylim(-2,2)

#    plt.subplot(133, aspect="equal")
#    plt.plot(x[2][indices], y[2][indices])
#    plt.scatter(0,0, marker="+")
#    plt.title("Not at resonance")
#    plt.xlabel(r"$x$ (along bar)")
#    plt.ylabel(r"$y$ (perpendicular to bar")
#    plt.xlim(-2,2)
#    plt.ylim(-2,2)
    
    plt.tight_layout()
    
    return

interact(plot_solution, plot_start=(0,plot_step_size*np.floor(t[0,-1] / plot_step_size-1),plot_step_size));


    

    




plot range: t=  0.0  -  49.0






    

In [5]:

    

plt.plot(x[0,:], y[0,:])
plt.title("Corotation Resonance")
plt.xlabel(r"$x$ (along bar)")
plt.ylabel(r"$y$ (perpendicular to bar")


    

    
Out[5]:





<matplotlib.text.Text at 0x7f7fb5420d90>






    

In [6]:

    

plt.plot(x[1,:], y[1,:])
plt.title("Outer Lindblad Resonance")
plt.xlabel(r"$x$ (along bar)")
plt.ylabel(r"$y$ (perpendicular to bar")


    

    
Out[6]:





<matplotlib.text.Text at 0x7f7fafc4a750>






    

In [7]:

    

plt.plot(x[2,:], y[2,:])
plt.title("No Resonance")
plt.xlabel(r"$x$ (along bar)")
plt.ylabel(r"$y$ (perpendicular to bar")


    

    
Out[7]:





<matplotlib.text.Text at 0x7f7fafc2c0d0>