第3章 重　力 　 【リスト3-1】重力ポテンシャルの表示 from tkinter import * import math,time #■初期設定 NumX = 50; NumY = 50 X0=250; Y0=250; Alfa=math.pi/9; Beta=math.pi/9 Z=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ID=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ZX=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ZY=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] #■隠れ線処理の初期化 def initHidden(): global X0, Y0, Alfa, Beta XY=[[0,0] for i in range(4)] XY[0][0]=400;XY[0][1]=400; XY[1][0]=400;XY[1][1]=410 XY[2][0]=410;XY[2][1]=410; XY[3][0]=410;XY[3][1]=400 DX=5 cosA=DX*math.cos(Alfa);cosB=DX*math.cos(Beta) sinA=DX*math.sin(Alfa);sinB=DX*math.sin(Beta) for i in range(NumX+1): for j in range(NumY+1): ZX[i][j]=X0+(i*cosA-j*cosB) ZY[i][j]=Y0-(i*sinA+j*sinB) for i in range(NumX): for j in range(NumY): ID[i][j]=canvas.create_polygon(XY,fill="white", outline="#000077") #■隠れ線処理 def hiddenLine(): #　画面の場合のみ可能（プロッタでは使えない） global canvas, NumX, NumY XY=[0 for i in range(8)] for i in range(NumX): i2=NumX-i;i1=i2-1 for j in range(NumY): if i>NumX/2 and j>NumY/2: continue #前方をカット j2=NumX-j; j1=j2-1 XY[0]=ZX[i1][j1]; XY[1]=ZY[i1][j1]-Z[i1][j1] XY[2]=ZX[i1][j2]; XY[3]=ZY[i1][j2]-Z[i1][j2] XY[4]=ZX[i2][j2]; XY[5]=ZY[i2][j2]-Z[i2][j2] XY[6]=ZX[i2][j1]; XY[7]=ZY[i2][j1]-Z[i2][j1] canvas.coords(ID[i][j],XY) Z=[[0 for j in range(NumY+1)] for i in range(NumX+1)] Center=NumX/2 # 質点の位置 r=2 # 質点の半径 rr=r*r # 質点の半径の２乗 c=200 # 表示上の重力定数 CR=-c/math.sqrt(rr) # 質点の表面の重力ポテンシャル collision= False #衝突の8定結果 def setData(): #データを設定 global canvas, NumX, NumY, Center, r,rr, c, CR for j in range(NumX+1): dj=j-Center; jj=dj*dj for k in range(NumY+1): dk=k-Center; kk=dk*dk; rjk=jj+kk if rjkrr: CCR =-c/math.sqrt(rjk) R[1]-=CCR return R def drawBall(X,Y):#重力場上の位置にボールを描く global canvas,BID R=screenXY(X,Y) BID=canvas.create_oval(R[0]-4,R[1]-4, R[0]+4,R[1]+4, fill="red") def drawBall2(X,Y):#Z=0の位置にボールを描く global canvas,BID2 R=screenXY2(X,Y) BID2=canvas.create_oval(R[0]-4,R[1]-4, R[0]+4,R[1]+4, fill="yellow") def moveBall(BX,BY):#ボール位置を変更(アニメーション用) global canvas,BID, BID2 R=screenXY(BX,BY); canvas.coords(BID, R[0]-4,R[1]-4, R[0]+4,R[1]+4) R=screenXY2(BX,BY); canvas.coords(BID2, R[0]-4,R[1]-4, R[0]+4,R[1]+4) #■メイン処理 root=Tk();root.title("重力ポテンシャル");root.geometry("500x300") canvas=Canvas(root,width=500,height=300,bg="white"); canvas.pack() initHidden(); setData();hiddenLine();root.update() root.title("重力ポテンシャル G= %7.4f, V0 = %7.4f" %(gFact, Vx)) drawBall(BX,BY);drawBall2(BX,BY) root.update() while(BX>=0 and BX<=NumX and BY>=0 and BY<=NumY): if collision :print("衝突発生"); break time.sleep(0.01) BX+=Vx; BY+=Vy drawBall(BX,BY);drawBall2(BX,BY) ACC=grad(BX,BY) Vx+=ACC[0]*gFact;Vy+=ACC[1]*gFact root.update() ３．３　質点別引力によるシミュレーション 【リスト3-3】質点別引力によるシミュレーション from tkinter import * import math, time X0=150; Y0=150;SC=1/5; RS=1/6 def genPlanet(canvas, X,Y,R,C):# 惑星表示生成 global X0,Y0,SC,RS; RR=RS*R; XX=X*SC+X0; YY=Y*SC+Y0 return canvas.create_oval(XX-RR, YY-RR,XX+RR, YY+RR, fill=C) def movePlanet(canvas, ID, X,Y,R):# 惑星位置変更 global X0,Y0,SC,RS; RR=RS*R; XX=X*SC+X0; YY=Y*SC+Y0 canvas.coords(ID, XX-RR, YY-RR,XX+RR, YY+RR) class Planet:#惑星クラス def __init__(self, canvas, R,M, C="red", X = 0, Y=0, DX=0, DY=0, AX=0,AY=0): self.canvas=canvas self.R = R # 表示半径 self.M = M # 質量 self.X = X # X座標 self.Y = Y # Y座標 self.DX = DX # X方向速度 self.DY = DY # Y方向速度 self.AX = AX # X方向加速度 self.AY = AY # Y方向加速度 self.C = C # 表示カラー self.ID = genPlanet(canvas, X,Y,R,C) Cfact=1; coll=False;DT=2 def calc():#計算 global numP,planet, Cfact,DT for i in range(numP): XX = planet[i].X; YY = planet[i].Y AX = 0; AY = 0 for j in range(numP): if i!=j: DX = planet[j].X - XX; DX2 = DX * DX DY = planet[j].Y - YY; DY2 = DY * DY M1 = planet[j].M * Cfact R2 = DX2 + DY2; R = math.sqrt(R2) CX = 1; CY=1 if DX < 0 : CX = -1 if DY < 0 : CY = -1 if DX2 < 0.000001:#衝突の判定 if DY2 < 0.000001 : coll=True; return AY = AY + M1 * CY / DY2 elif DY2 < 0.000001: AX = AX + M1 * CX / DX2 else: AA = M1 / R2 AX = AX + AA * DX / R; AY = AY + AA * DY / R planet[i].AX = AX; planet[i].AY = AY planet[i].DX = planet[i].DX + AX * DT planet[i].DY = planet[i].DY + AY * DT planet[i].X = planet[i].X + planet[i].DX * DT planet[i].Y = planet[i].Y + planet[i].DY * DT pile=False # 惑星を上描き表示するかどうかのフラグ def dsp():#表示 global numP,planet,canvas,heapUp for i in range(numP): if pile: genPlanet(canvas, planet[i].X, planet[i].Y, planet[i].R, planet[i].C) else: movePlanet(canvas, planet[i].ID, planet[i].X, planet[i].Y, planet[i].R) #■メイン処理 root=Tk();root.title("惑星運動シミュレーション") canvas=Canvas(root, width=300,height=300,bg="white") canvas.pack() numP=4 planet=[0 for i in range(numP)] planet[0]=Planet(canvas,50,10000,"red") planet[1]=Planet(canvas,30,10,"skyblue", 600, 0, 0, 4) planet[2]=Planet(canvas,15, 5,"yellow" , 0, 250, -6, 0) planet[3]=Planet(canvas,20, 5,"green" , 0, 400, -5, 0) if pile : n=0 while(not coll): if pile: n+=1 #上描き表示のとき固定回数とする if n>100: break root.update() time.sleep(0.01) calc() dsp() 【初期設定の変更】 　 numP=4 planet=[0 for i in range(numP)] planet[0]=Planet(canvas,50,10000,"red") planet[1]=Planet(canvas,30,10,"skyblue" , 600, 0, 0, 4) planet[2]=Planet(canvas,15, 0.1,"yellow", 640, 0, -2, 2) planet[3]=Planet(canvas,20, 5,"green" , 0, 400, -5, 0) 　 【2体だけの例】 Cfact=1; coll=False;DT=4 numP=2 planet=[0 for i in range(numP)] planet[0]=Planet(canvas,20,1000,"red" , 40, 0, 2.2, 0.8) planet[1]=Planet(canvas,20,1000,"yellow", 0, 160, -2.2, -0.8)