66 lines
No EOL
1.4 KiB
Python
66 lines
No EOL
1.4 KiB
Python
#pylint: disable=W0401,W0614,W0622
|
|
#All pylint codes: http://pylint.pycqa.org/en/latest/technical_reference/features.html
|
|
|
|
import numpy as np
|
|
import ctypes, subprocess
|
|
|
|
recompile = False
|
|
|
|
if recompile:
|
|
p = subprocess.Popen('g++ -O3 -I/home/meiron/boost_1_72_0 -shared -o libmain.so -fPIC loadtxt.cpp main.cpp -lgsl'.split())
|
|
p.wait()
|
|
if p.returncode != 0: raise RuntimeError(p.returncode)
|
|
|
|
libmain = ctypes.CDLL('./libmain.so')
|
|
def integrate(y0, t_max, step_size=None):
|
|
y0 = (ctypes.c_double*6)(*y0)
|
|
if step_size is None: step_size = t_max
|
|
size = int(t_max // step_size) + 1
|
|
y = (ctypes.c_double*size*6)()
|
|
libmain.integrate(y0, ctypes.c_double(t_max), ctypes.c_double(step_size), y)
|
|
y = np.array(y).reshape(size, 6)
|
|
return np.array(y)
|
|
|
|
from pylab import *
|
|
t_max = 10
|
|
x_max = 120
|
|
ic = array([120,0,0,0,30,0])
|
|
|
|
step_size = 32/4096
|
|
res = integrate(ic, t_max, step_size=step_size)
|
|
|
|
x, y, z, vx, vy, vz = res.T
|
|
t = arange(0,t_max+step_size,step_size)
|
|
|
|
figure(figsize=(16,16))
|
|
subplot(221)
|
|
plot(x, y)
|
|
plot(x[0], y[0], 'x')
|
|
gca().set_aspect('equal')
|
|
xlabel('x')
|
|
ylabel('y')
|
|
xlim(-x_max,x_max)
|
|
ylim(-x_max,x_max)
|
|
|
|
subplot(223)
|
|
plot(x, z)
|
|
plot(x[0], z[0], 'x')
|
|
gca().set_aspect('equal')
|
|
xlabel('x')
|
|
ylabel('z')
|
|
xlim(-x_max,x_max)
|
|
ylim(-x_max,x_max)
|
|
|
|
subplot(224)
|
|
plot(y, z)
|
|
plot(y[0], z[0], 'x')
|
|
gca().set_aspect('equal')
|
|
xlabel('z')
|
|
ylabel('y')
|
|
xlim(-x_max,x_max)
|
|
ylim(-x_max,x_max)
|
|
|
|
subplots_adjust(wspace=0.25)
|
|
|
|
savefig('test_orbit.png')
|
|
show() |