python勉強中

なんだか新しいプログラム言語をやってみたくなって、pythonを勉強中。
以前にjavascriptで書いた太陽系形成プログラムを書き直してみた。行列演算が出来るのでforループを書く必要が殆どない。むしろループを回さずにnativeの行列演算ルーチンを使うとMATLAB同様とても速い。MATLABよりアニメーションかはやりやすいようだ。MP4ではき出すのが簡単にできる。フリーなのも有り難い。説明書がいろんな所に散らばっていて、読みにくい感じがするのがちょっと気になる。

それから粒子法による物理シミュレーションのやり方も勉強中なので、圧力とか温度とかも入れて中心の塊を核融合で輝かせてみたい。

www.youtube.com

import numpy as np
import random
import copy
from matplotlib import pyplot as plt
import matplotlib.animation as animation
fig = plt.figure()
ims=[]

MassTotal=100
N=100
dt=0.1
h=0.3
G=0.0001*MassTotal/N
tau=20.0
x=(np.random.rand(N)-0.5)*10
y=(np.random.rand(N)-0.5)*10
X=np.random.rand(N)
Y=np.random.rand(N)
v=np.sqrt(-2.0*np.log(X))*np.cos(2*np.pi*Y)
theta=np.pi*np.random.rand(N)
vx=v*np.cos(theta)/20
vy=v*np.sin(theta)/20
for i in range(20000):
    dx=x[np.newaxis,:].T-x
    dy=y[np.newaxis,:].T-y
    dvx=vx[np.newaxis,:].T-vx
    dvy=vy[np.newaxis,:].T-vy
    r2=dx*dx+dy*dy
    r=np.sqrt(r2)
    r3=r2*r+h**4
    temp=r<h
    f0=G*(~temp*1.0/r3+temp*1.0)
    fx=np.sum(f0*dx,axis=0)+np.sum(dvx*temp/tau,axis=0)
    fy=np.sum(f0*dy,axis=0)+np.sum(dvy*temp/tau,axis=0)
    x+=vx*dt
    y+=vy*dt
    vx+=fx*dt
    vy+=fy*dt
    if i%10==0:
        frame,=plt.plot(x,y,"bo",markersize=2)
        plt.xlim((-20,20))
        plt.ylim((-20,20))
        plt.gca().set_aspect('equal', adjustable='box')
        ims.append([frame])

ani = animation.ArtistAnimation(fig, ims, interval=5, repeat_delay=1000)
ani.save('im.mp4', writer="ffmpeg",fps=30)
plt.show()