带电粒子在电磁场中的运动的计算机模拟
样例
样例1
模型描述：在一个布满整个空间、垂直纸面向下的匀强磁场中，两个等质量、带等量电荷但电荷量相反的两个带电粒子由静止释放。
指令:
md1.txt:
adm
-7500
True
N
-10000000
-10000000
10000000
-10000000
10000000
10000000
-10000000
10000000
end
adb
9.5
-9
1
0
0
0.0003
adb
9.5
9
1
0
0
-0.0003
s
10
0.01
1
Y
然后输入python SMCPEF.py md1 <md1.txt
样例2
模型描述：一个带电粒子在一个布满整个空间、沿纸面向下的匀强电场中作类斜抛运动。
样例3
模型描述：一个带电粒子在一个最普通的跑道形回旋加速器中的运动
指令:
md3.txt
adb
0
-9
1
0
0
1
ade
3
0
-1<x<1 and -9.1<y<-8.9
N
-1
-9.1
1
-9.1
1
-8.9
-1
-8.9
end
adm
-1
x>=1
N
1
-100000
100000
-100000
100000
100000
1
100000
end
adm
-1
x<=-1
N
-1
100000
-100000
100000
-100000
-100000
-1
-100000
end
s
10
0.001
1
Y
然后输入python SMCPEF.py md3 <md3.txt
源代码
SMCPEF.py
import os
import sys
import matplotlib.pyplot as plt
import time
import math
import numpy as np
import cv2
from math import sin,cos,tan,pi,sqrt
if len(sys.argv)==1:
while True:
s=input("请输入项目名称(不能有空格) >>>")
os.system("python "+sys.argv[0]+" "+s)
MyName=sys.argv[1]
EBalls=[]
FMags=[]
FEles=[]
Boards=[]
Reflcs=[]
Geners=[]
ori_pics=[]
CameraMove=False
Lx=-10
Ly=10
Width=20
#----------------------------------------------
e_k=8987551787
TC=0
plt.figure(figsize=(10.24,10.24))
#----------------------------------------------
class Pic_t:
def __init__(self,ti):
self.time=ti
self.obj=[]
def add_object(self,o):
self.obj.append(o)
class Vector:
def __init__(self,x,y):
self.x=x
self.y=y
def Add(A,B):
return Vector(A.x+B.x,A.y+B.y)
def Minus(A,B):
return Vector(A.x-B.x,A.y-B.y)
def NMt(k,A):
return Vector(A.x*k,A.y*k)
class Ball:
def __init__(self,name,X,Y,m,vx,vy,q):
self.name=name
self.X=X
self.Y=Y
self.m=m
self.vx=vx
self.vy=vy
self.q=q
self.alive=1
def Draw_Bound(self):
pass
def Force(self,t,bb):
global e_k
A=Vector(self.X,self.Y)
B=Vector(bb.X,bb.Y)
t=-e_k*self.q*bb.q
C=Minus(A,B)
FS=t/(C.x*C.x+C.y*C.y)
L=math.sqrt(C.x*C.x+C.y*C.y)
C=NMt(FS/L,C)
return C
class Mag:
def __init__(self,name):
self.name=name
def Draw_Bound(self):
pass
def B(self,t,x,y):
pass
def inm(self,x,y):
pass
def Force(self,t,ball):# F=qvB
if not self.inm(ball.X,ball.Y):
return Vector(0,0)
b=self.B(t,ball.X,ball.Y) # B=(0,0,b) v=(vx,vy,0) (vy*b,-vx*b,0)
return NMt(ball.q,Vector(ball.vy*b,-ball.vx*b))
class Ele():
def __init__(self,name):
self.name=name
def Draw_Bound(self):
pass
def E(self,t,x,y):
pass
def ine(self,x,y):
pass
def Force(self,t,ball):
if not self.ine(ball.X,ball.Y):
return Vector(0,0)
v=self.E(t,ball.X,ball.Y)
v=NMt(ball.q,v)
return v
class Board:
def __init__(self,name,a,b,c,d):
self.name=name
self.va=a
self.vb=b
self.vc=c
self.lm=d
def Draw_Bound(self):
pass
def itr(self,x1,y1,x2,y2):
pass
def idv(self,x,y):
pass
def bounce(self,ball):
pass
class LBoard(Board):
def __init__(self,name,A,B,C,L):
super(Board,self).__init__(name,A,B,C,L)
def Draw_Bound(self):
pass
def itr(self,x1,y1,x2,y2):
pass
def idv(self,x,y):
pass
def bounce(self,ball):
pass
class CBoard(Board):
def __init__(self,name,x,y,R,L):
super(Board,self).__init__(name,x,y,R,L)
def Draw_Bound(self):
pass
def itr(self,x1,y1,x2,y2):
pass
def idv(self,x,y):
pass
def bounce(self,ball):
pass
def Engine():# 物理引擎
global TC,TR
#print(len(FMags))
while TC<=TR:
#--------------------------
dt=Tik(TC)
totBall=len(EBalls)
for i in range(totBall):
dx=NMt(dt,Vector(EBalls[i].vx,EBalls[i].vy))
#---------------------------判断碰撞
#---------------------------
EBalls[i].X+=dx.x
EBalls[i].Y+=dx.y
F=Vector(0,0)
for j in range(totBall):
if i!=j:
F=Add(F,EBalls[j].Force(TC,EBalls[i]))
for j in FMags:
F=Add(F,j.Force(TC,EBalls[i]))
#print(F.x,F.y)
for j in FEles:
F=Add(F,j.Force(TC,EBalls[i]))
dv=NMt(dt/EBalls[i].m,F)
EBalls[i].vx+=dv.x
EBalls[i].vy+=dv.y
#print(EBalls[i].name,EBalls[i].X,EBalls[i].Y)
#--------------------------
nP=Pic_t(TC)
for i in EBalls:
nP.add_object(["Ball",i.name,i.X,i.Y])
#print("length",len(nP.obj))
ori_pics.append(nP)
TC=TC+dt
def NeedMove(point):
global Lx,Ly,Width
tot=len(ori_pics[point].obj)
for i in range(tot):
if ori_pics[point].obj[i][0]=="Ball":
if Lx<=ori_pics[point].obj[i][2]<=Lx+Width and Ly-Width<=ori_pics[point].obj[i][3]<=Ly:
pass
else:
#print("2!!!",item.X,item.Y)
return 2
for i in range(tot):
if ori_pics[point].obj[i][0]=="Ball":
for j in range(tot):
if i!=j and ori_pics[point].obj[j][0]=="Ball":
dx=ori_pics[point].obj[i][2]-ori_pics[point].obj[j][2]
dy=ori_pics[point].obj[i][3]-ori_pics[point].obj[j][3]
if sqrt(dx*dx+dy*dy)<=0.01*Width:
print("1!!!")
return 1
return 0
def PicBuilder():#20帧的视频
global BS,ori_pics,Lx,Ly,Width
realtimepertik=0.05*BS
point=0
tikid=0
while point<len(ori_pics):
Rtime=realtimepertik*tikid
while point<len(ori_pics):
if point+1<len(ori_pics) and ori_pics[point+1].time<=Rtime:
point+=1
else:
break
if point==len(ori_pics)-1:
break
#开始画
if CameraMove==True:
st=NeedMove(point)
if st!=0:
Minx=10.0**20.0
Maxx=-Minx
Miny=Minx
Maxy=-Minx
for item in ori_pics[point].obj:
if item[0]!="Ball":
continue
Minx=min(Minx,item[2])
Maxx=max(Maxx,item[2])
Miny=min(Miny,item[3])
Maxy=max(Maxy,item[3])
ttt=max(Maxx-Minx,Maxy-Miny)
if ttt==0:
ttt=19/2
Width=int(2*ttt+1)
Lx=Minx-0.25*Width
Ly=Maxy+0.25*Width
else:
pass
#print([Lx,Lx+Width,Ly-Width,Ly])
plt.axis([Lx,Lx+Width,Ly-Width,Ly])
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
#plt.scatter(x_values, y_values, s=100)
#print("size:",len(ori_pics[point].obj))
for item in FMags:
item.Draw_Bound()
for item in FEles:
item.Draw_Bound()
for item in ori_pics[point].obj:
if item[0]=="Ball":
plt.scatter(item[2],item[3],color="black")
#print(item[2],item[3])
plt.savefig("tmp/"+str(tikid)+".png")
plt.clf()
tikid=tikid+1
return tikid-1
def VideoBuilder(tot):
size = (1024,1024)
videowrite = cv2.VideoWriter(MyName+'.mp4',cv2.VideoWriter_fourcc('M', 'P', '4', 'V'),20,size)
img_array=[]
for filename in [r'tmp\\\{0}.png'.format(i) for i in range(tot+1)]:
img = cv2.imread(filename)
if img is None:
print(filename + " is error!")
continue
img_array.append(img)
for i in range(tot+1):
videowrite.write(img_array[i])
videowrite.release()
def Render():
a=PicBuilder()
print("照片生成完毕，正在合成视频")
VideoBuilder(a)
#--------------------------------------------------------
TR,BS,Tik_time=0,0,0
block_len=20
def Tik(T):
return Tik_time
def G1(T):
gtime=2
if T>=gtime:
pass
Geners.append(G1)
Geners[0](TC)
def Quit():
quit(0)
def Adv():
os.system("copy "+sys.argv[0]+" 高级版本.py")
input("请按打开 高级版本.py 自行编辑 按Enter退出当前会话")
def Adb():
global EBalls
#def __init__(self,name,X,Y,m,vx,vy,q):
name=""#input("请输入小球名称")
X=float(input("请输入小球X坐标"))
Y=float(input("请输入小球Y坐标"))
m=float(input("请输入小球质量"))
vx=float(input("请输入小球X方向分速度"))
vy=float(input("请输入小球Y方向分速度"))
q=float(input("请输入小球电荷量"))
EBalls.append(Ball(name,X,Y,m,vx,vy,q))
class EF1(Ele):
def __init__(self,name,Ex,Ey,inf,X,Y):
Ele.__init__(self,name)
self.name=name
self.Ex=Ex
self.Ey=Ey
self.inf=inf
self.X=X
self.Y=Y
def Draw_Bound(self):
plt.fill(self.X,self.Y, 'r')
def E(self,t,x,y):
t=t
x=x
y=y
ex=eval(self.Ex)
ey=eval(self.Ey)
return Vector(ex,ey)
def ine(self,x,y):
x=x
y=y
return eval(self.inf)
def Ade():
name=""#input("请输入电场名称")
Ex=input("请输入点(x,y)处在t时刻场强的X分量（用无赋值的单行python表达式）")
Ey=input("请输入点(x,y)处在t时刻场强的Y分量（用无赋值的单行python表达式）")
inf=input("请输入点(x,y)是否位于电场内部的判别式")
s=input("是否为圆形? Y/N ")
X=[]
Y=[]
if s=="Y":
x=float(input("X"))
y=float(input("Y"))
r=float(input("R"))
for i in range(1000):
theta=2*pi*i/1000
X.append(x+r*cos(theta))
Y.append(y+r*sin(theta))
else:
print("接下来请依次逆时针输入边界坐标，无需重复输入第一个点，结束时输入end")
tot=0
while True:
tot+=1
s=input("第"+str(tot)+"个的X坐标")
if s=="end":
break
X.append(float(s))
s=input("第"+str(tot)+"个的Y坐标")
if s=="end":
break
Y.append(float(s))
FEles.append(EF1(name,Ex,Ey,inf,X,Y))
class MF1(Mag):
def __init__(self,name,Bs,inf,X,Y):
Mag.__init__(self,name)
self.name=name
self.Bs=Bs
self.inf=inf
self.X=X
self.Y=Y
def Draw_Bound(self):
plt.fill(self.X,self.Y, 'blue')
def B(self,t,x,y):
return eval(self.Bs)
def inm(self,x,y):
return eval(self.inf)
def Adm():
name=""#input("请输入磁场名称")
B=input("请输入点(x,y)处在t时刻场强(正为指向纸外)（用无赋值的单行python表达式）")
inf=input("请输入点(x,y)是否位于磁场内部的判别式")
s=input("是否为圆形? Y/N ")
X=[]
Y=[]
if s=="Y":
x=float(input("X"))
y=float(input("Y"))
r=float(input("R"))
for i in range(1000):
theta=2*pi*i/1000
X.append(x+r*cos(theta))
Y.append(y+r*sin(theta))
else:
print("接下来请依次逆时针输入边界坐标，无需重复输入第一个点，结束时输入end")
tot=0
while True:
tot+=1
s=input("第"+str(tot)+"个的X坐标")
if s=="end":
break
X.append(float(s))
s=input("第"+str(tot)+"个的Y坐标")
if s=="end":
break
Y.append(float(s))
FMags.append(MF1(name,B,inf,X,Y))
def Start():
global TR,Tik_time,BS,CameraMove,Lx,Ly,Width
TR=float(input("渲染总时长"))
Tik_time=float(input("单步时间间隔"))
BS=float(input("视频倍速"))
cm=input("是否打开摄像头自动移动？(否则将为恒定范围）Y/N")
if cm=='Y':
CameraMove=True
else:
Lx=float(input("请输入左上角X坐标"))
Ly=float(input("请输入左上角Y坐标"))
Width=float(input("请输入视野宽度"))
Engine()
print("物理运算完毕，开始渲染")
print("长度：",len(ori_pics))
Render()
input("合成完成，请按Enter推出并查看"+MyName+".mp4")
Quit()
Menu={"adv":Adv,"s":Start,"adb":Adb,"ade":Ade,"adm":Adm,"q":Quit}
while True:
s=input("""请输入命令：
adv 进入高级模式
s 开始模拟
adb 增加一个带电小球
ade 增加一个电场
adm 增加一个磁场
q 退出
>>> """)
if s in ["adv","s","adb","ade","adm","q"]:
Menu[s]()
本作品由happyZYM采用知识共享 署名-非商业性使用-相同方式共享 4.0 (CC BY-NC-SA 4.0) 国际许可协议(镜像(简单版)镜像(完整版))进行许可。
转载请注明出处：https://www.cnblogs.com/happyZYM/p/SMCPEF.html (近乎)全文转载而非引用的请在文首添加出处链接。
