作成した3次元グラフの座標を出力させたい

Question

実現したいこと

matplotlibで3次元のグラフを作成しました。
ソースコードでは、４つの点の軌道が描かれており、４つの点が時刻tの時のx座標、y座標、z座標をprintで出力させたいです。
時刻tに関しては始点から終点までの時間であり、始点から終点までのx座標、y座標、z座標をprintで出力させる方法を教えていただきたいです。
例えば、
r1_sol
t=1, x=0, y=0, z=0
t=2, x=1, y=1, z=1
...
r2_sol
t=1, x=0, y=0, z=0
t=2, x=1, y=1, z=1
...
r3_sol
t=1, x=0, y=0, z=0
t=2, x=1, y=1, z=1
...
r4_sol
t=1, x=0, y=0, z=0
t=2, x=1, y=1, z=1
...
のような感じです。

該当のソースコード

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

#万有引力定数を定義
G=6.67e-11 #N-m2/kg2
#参照数量
m_nd=1.989e+30 #kg #太陽の質量
r_nd=5.326e+12 #m #アルファケンタウリの星間距離
v_nd=30000 #m/s #太陽の周りの地球の相対速度
t_nd=79.91*365*24*3600*0.51 #s #アルファケンタウリの軌道期間
#正味定数
K1=G*t_nd*m_nd/(r_nd**2*v_nd)
K2=v_nd*t_nd/r_nd

#質量の定義
m1=1 #質量が太陽の５倍の質点
m2=1 #質量が太陽の質点
#位置の初期条件
r1=[-0.5,0,0] #m
r2=[0.5,0,0] #m
r1=np.array(r1,dtype="float64")
r2=np.array(r2,dtype="float64")
r_com=(m1*r1+m2*r2)/(m1+m2)
#速度の初期条件
v1=[0.01,0.01,0] #m/s
v2=[-0.05,0,-0.1] #m/s
v1=np.array(v1,dtype="float64")
v2=np.array(v2,dtype="float64")
v_com=(m1*v1+m2*v2)/(m1+m2)

#運動方程式を定義
def TwoBodyEquations(w,t,G,m1,m2):
    r1=w[:3]
    r2=w[3:6]
    v1=w[6:9]
    v2=w[9:12]
    r=np.linalg.norm(r2-r1)
    dv1bydt=K1*m2*(r2-r1)/r**3
    dv2bydt=K1*m1*(r1-r2)/r**3
    dr1bydt=K2*v1
    dr2bydt=K2*v2
    r_derivs=np.concatenate((dr1bydt,dr2bydt))
    derivs=np.concatenate((r_derivs,dv1bydt,dv2bydt))
    return derivs

#パッケージの初期パラメータ
init_params=np.array([r1,r2,v1,v2]) #create array of initial params
init_params=init_params.flatten() #flatten array to make it 1D
time_span=np.linspace(0,8,500) #8 orbital periods and 500 points
#odeを実行
from scipy.integrate import odeint
two_body_sol=odeint(TwoBodyEquations,init_params,time_span,args=(G,m1,m2))

r1_sol=two_body_sol[:,:3]
r2_sol=two_body_sol[:,3:6]

rcom_sol=(m1*r1_sol+m2*r2_sol)/(m1+m2)
r1com_sol=r1_sol-rcom_sol
r2com_sol=r2_sol-rcom_sol

#図を作成
fig=plt.figure(figsize=(5,5))
ax=fig.add_subplot(111,projection="3d")
ax.plot(r1_sol[:,0],r1_sol[:,1],r1_sol[:,2],color="blue")
ax.plot(r2_sol[:,0],r2_sol[:,1],r2_sol[:,2],color="red")
ax.plot(r3_sol[:,0],r3_sol[:,1],r3_sol[:,2],color="green")
ax.plot(r4_sol[:,0],r4_sol[:,1],r4_sol[:,2],color="orange")
plt.savefig("FourBodyEquation10.pdf")
#さらにいくつかのベルとホイッスルを追加する
ax.set_xlabel("x-coordinate",fontsize=5)
ax.set_ylabel("y-coordinate",fontsize=5)
ax.set_zlabel("z-coordinate",fontsize=5)

steps = r1_sol.shape[0]
artists = [
    [ax.scatter(*np.c_[r1_sol[i], r2_sol[i], r3_sol[i], r4_sol[i]],
                color=["blue", "red","green","orange"], marker="o", s=30)]
    for i in range(0, steps, 2)
]
ani = animation.ArtistAnimation(fig, artists, interval=50)
ani.save('FourBodyEquations2.gif', writer='imagemagick')``

Comments

[y_waiwai]

時間の関数としての座標、ってだけでは意味不明なので、もうちょい具体的に説明しよう
実例も出していただければ。

[jbpb0]

質問のコードをそのまま実行したら、エラーになります

NameError: name 'r3_sol' is not defined

NameError: name 'r4_sol' is not defined

ちゃんと動くのを掲載した方が、他人が現象を再現しやすくなるので、回答されやすくなると思いますよ

Answer 1

name = 'r1_sol'
exec("rtmp_sol = %s" % name)
print(name)
for i in range(rtmp_sol.shape[0]):
    print('t=' + str(i+1) + ', x=' + str(rtmp_sol[i, 0]) + ', y=' + str(rtmp_sol[i, 1]) + ', z=' + str(rtmp_sol[i, 2]))

「name =」の後を「'r2_sol'」とかに変えて実行してください

Comments

[osamuhito]

ありがとうございます！