%display latex #latexを使って、結果を表示します。
load( "setupSphericalCoordinates.sage")
E.<x,y,z> = EuclideanSpace(name="E")
show(E.atlas(),E.frames(),[type(f) for f in E.frames()])
SPC, SPF=E.spherical_coordinates(), E.spherical_frame(); # E.spherical_frame()は 正規直交基底を与える。こ
CC, CF= E.cartesian_coordinates(), E.cartesian_frame()
NF=E.frames()[0]
# metric tensor
g=E.metric()
# covariant derivative operator.
nab=g.connection(init_coef=True) #init_coef=True is a default value
E.default_frame(),E.default_chart()
g[NF,:], g.up(g).up(g)[NF,:]
CC.frame(),E.frames()
E.changes_of_frame()
E.change_of_frame(SPF,CF)[:],E.change_of_frame(CF,SPF)[:],
import numpy as np
import matplotlib.pyplot as pyplot
w = 1
Z,X = np.mgrid[-w*64:w*64+1:1, -w*64:w*64+1:1]/64
gamma=2
R= np.sqrt(gamma**2*X**2+Z**2)
Q=0.001
U = X/R
V = Z/R
speed = np.sqrt(U**2 + V**2)
fig, axs = pyplot.subplots(3, 2, figsize=(7, 9), height_ratios=[1, 1, 2])
axs = list(axs.flat)
# Varying density along a streamline
axs[0].streamplot(X, Z, U, V, density=[0.5, 1])
axs[0].set_title('Varying Density')
# Varying color along a streamline
strm = axs[1].streamplot(X, Z, U, V, color=U, linewidth=2, cmap='autumn')
fig.colorbar(strm.lines)
axs[1].set_title('Varying Color')
# Varying line width along a streamline
lw = np.sqrt(np.sqrt(speed))
axs[2].streamplot(X, Z, U, V, density=0.6, color='k', linewidth=lw, broken_streamlines=True)
axs[2].set_title('Varying Line Width')
axs[5].streamplot(X, Z, U, V, broken_streamlines=False)
axs[5].set_title('Streamplot with unbroken streamlines')
# Create a mask
mask = np.zeros(U.shape, dtype=bool)
#mask[40:60, 40:60] = True
U[:20, :20] = np.nan
U = np.ma.array(U, mask=mask)
axs[4].streamplot(X, Z, U, V, color='r')
axs[4].set_title('Streamplot with Masking')
pyplot.tight_layout()
pyplot.show()
var('gamma')
R=sqrt(gamma**2*x**2+(y+1e-6)**2+z**2)
EMQ=E.vector_field((gamma*x/R**3, gamma*(y+1e-6)/R**3, gamma*z/R**3),
frame=CF,chart=CC)
goptions=dict(chart=CC, frame=CF,
coords=(x,y,z),
ambient_coords=(x,y,z),
ranges= {x:(-1,1),y:(-1,1),z:(-1,1)},
number_values={x:4,y:4,z:8},
parameters={gamma:N(1/sqrt(1-0.8**2))},
width=1, scale=.6,)
cmap=colormaps["rainbow"]
rcmap=cmap.reversed()
plt=sum((
CC.plot(number_values={x:3,y:3,z:3}, chart=CC, frame=CF,
ambient_coords=(x,y,z),label_axes=False,
color='gray',linewidth=0.05, aspect_ratio=1,
ranges= {x:(-1,1),y:(-1,1),z:(-1,1)}),
EMQ.plot( color=rcmap(255//3) , label_axes=False, **goptions) ,
#EMQ.plot( color=rcmap(2*255//3) , label_axes=False, **goptions),
#EMQ.plot( color=rcmap(255), label_axes=False, **goptions),
))
plt.show(frame=False)
var('gamma')
R=sqrt(gamma**2*x**2+y**2+z**2+1e-6)
EMQ=E.vector_field((x/R*log(R**2), log(R**2)*y/R**3, log(R**2)*z/R**3),
frame=CF,chart=CC)
goptions=dict(chart=CC, frame=CF,
coords=(x,y,z),
ambient_coords=(x,y,z),
ranges= {x:(-1,1),y:(-1,1),z:(-1,1)},
number_values={x:4,y:4,z:8},
parameters={gamma:N(1/sqrt(1-0.8**2))},
width=1, scale=.6,)
cmap=colormaps["rainbow"]
rcmap=cmap.reversed()
plt=sum((
CC.plot(number_values={x:3,y:3,z:3}, chart=CC, frame=CF,
ambient_coords=(x,y,z),label_axes=False,
color='gray',linewidth=0.05, aspect_ratio=1,
ranges= {x:(-1,1),y:(-1,1),z:(-1,1)}),
EMQ.plot( color=rcmap(255//3) , label_axes=False, **goptions) ,
#EMQ.plot( color=rcmap(2*255//3) , label_axes=False, **goptions),
#EMQ.plot( color=rcmap(255), label_axes=False, **goptions),
))
plt.show(frame=False)
異なる$\gamma$での電場強度の等高線を描いてみましょう。 速度が0となる $\gamma=1$ではもちろん球対称ですが、 光速に近い $\gamma=5/3$ では下の図のように、電場強度の分布は進行方向に対して、扁平に圧縮された分布と なっていることがわかります。
N_sample=200
gamma=1
plt=contour_plot(log(
max_symbolic(1e-6, min_symbolic(1e6, gamma*sqrt(x**2+y**2)/sqrt(gamma**2*x**2+y**2)**3)
)),
xrange=(-2,2),yrange=(-1,1),
plot_points=N_sample,fill=True,
cmap="jet",
)
plt.save("./EF_StaticQ.png")
show(plt)
#E_r
gamma=5/3
def EM(x,y):
return log(gamma*sqrt(x**2+y**2)/sqrt(gamma**2*x**2+(y+0.0001)**2)**3
)
plt=contour_plot(EM(x,y),
xrange=(-2,2),yrange=(-1,1),
plot_points=N_sample,fill=True,
cmap="jet",
)
plt.save("./EF_movingQ.png")
show(plt)
#E_r
gamma=1/sqrt(1-0.95**2)
def EM(x,y):
return log(gamma*sqrt(x**2+y**2)/sqrt(gamma**2*x**2+(y+0.0001)**2)**3
)
plt=contour_plot(EM(x,y),
xrange=(-2,2),yrange=(-1,1),
plot_points=N_sample,fill=True,
cmap="jet",
)
show(plt)
#E_x
gamma=1/sqrt(1-0.9**2)
plt=contour_plot(
sign(x)*log( sqrt(x**2)/sqrt(gamma**2*x**2+y**2)**3 ) ,
xrange=(-0.2,0.2),yrange=(-0.1,0.1),
plot_points=N_sample,fill=True,
cmap="jet_r",
)
show(plt)
#Ey
plt=contour_plot(
sign(y)*log(
sqrt(y**2)/sqrt(gamma**2*x**2+y**2)**3
),
xrange=(-0.2,0.2),yrange=(-0.1,0.1),
plot_points=N_sample,fill=True,
cmap="jet_r",
)
plt.save("./EFY_movingQ.png")
show(plt)
gamma=5/3
v= 4/5
plts=[]
for i in range(0,50):
t=0.2*(i+0.5)-3
R1=sqrt(gamma**2*(x-v*t)**2+(y+10^(-7))**2)
R2=sqrt(gamma**2*(x+v*t)**2+(y+10^(-7))**2)
Ex= (x-v*t)/R1**3 - (x+v*t)/R2**3
Ey= (y+10^(-7))*(1/R1**3 - 1/R2**3)
Bz= v*(y+10^(-7))*(1/R1**3 + 1/R2**3)
plt=contour_plot(
max_symbolic(-3,min_symbolic(3, log(Ex**2+Ey**2))),
xrange=(-6,6),yrange=(-1,1),
plot_points=N_sample,fill=True,
cmap="jet",
contours=[
-4.278607150491874, -1.290575941139906, 1.697455268212062,
4.685486477564032, 7.673517686916002, 10.66154889626797,
13.64958010561994, 16.63761131497191]
)
plt2=contour_plot(
sign(Bz)*max_symbolic(-3 ,min_symbolic(3,log(Bz**2))),
xrange=(-6,6), yrange=(-1,1),
plot_points=N_sample,fill=True,
cmap="jet_r",
contours=[
-17.90030162, -12.78592973, -7.67155784,
-2.55718595, 2.55718595, 7.67155784,
12.78592973, 17.90030162],
)
if ( mod(i,10) == 1):
plt.save(f"./movie/EF_movingCharges_frame{i:02d}.png")
plt2.save(f"./movie/MF_movingCharges_frame{i:02d}.png")
(plt2+plt).save(f"./movie/movingCharges_frame{i:02d}.png")
plts.append(plt2+plt)
anim=animate(plts)
anim.show(delay=50, iterations=5)
anim.save("./movie/movingCharges.mp4",show_path=True)
Animation saved to file /Users/noboru/src/JupyterNotebooks/SageMath/SageManifolds/movie/movingCharges.mp4.
anim.show(delay=10, iterations=5)
v= 0.005
gamma=1/sqrt(1-v**2)
plts=[]
dt=0.3
#dt=1.0
n=0.6/(v*dt)
for i in range(0,int(n)):
t=dt*(i+0.5) - 0.2/v
R1=sqrt(gamma**2*(x - v*t)**2+(y + 10^(-6))**2)
R2=sqrt(gamma**2*(x + v*t)**2+(y + 10^(-6))**2)
Ex= (x - v*t)/R1**3 - (x + v*t)/R2**3
Ey= (y+10^(-6))*(1/R1**3 - 1/R2**3)
Bz= v*(y+10^(-6))*(1/R1**3 + 1/R2**3)
plt=contour_plot(
log(Ex**2+Ey**2),
xrange=(-1.5, 1.5),
yrange=(-0.5, 0.5),
plot_points=N_sample,fill=True,
contours=[-1.67973392, 1.31869241, 4.31711875,
7.31554508 , 10.31397142, 13.31239775,
16.31082409, 19.30925042],
cmap="jet",
)
plt2=contour_plot(
sign(Bz)*log(Bz**2),
xrange=(-1.5, 1.5),yrange=(-0.5, 0.5),
plot_points=N_sample,fill=True,
contours=[
-9.65789873, -6.8984991, -4.13909946,
-1.37969982, 1.37969982, 4.13909946,
6.8984991, 9.65789873],
cmap="jet_r",
)
if (mod(i,30)==0):
(plt2+plt).save(f"./movie/SlowMovingCharges_frame{i:03d}.png")
plts.append(plt2+plt)
anim=animate(plts)
anim.save("./movie/SlowMovingCharges.mp4",show_path=True)
Animation saved to file /Users/noboru/src/JupyterNotebooks/SageMath/SageManifolds/movie/SlowMovingCharges.mp4.
anim.show(delay=10, iterations=5)
v= 0.005
gamma=1/sqrt(1-v**2)
t=10
R1=sqrt(gamma**2*(x - v*t)**2+(y + 10^(-6))**2)
R2=sqrt(gamma**2*(x + v*t)**2+(y + 10^(-6))**2)
Ex= (x - v*t)/R1**3 - (x + v*t)/R2**3
Ey= (y+10^(-6))*(1/R1**3 - 1/R2**3)
Bz= v*(y+10^(-6))*(1/R1**3 + 1/R2**3)
for range in (2,20,200,2000):
plt=contour_plot(
log(Ex**2+Ey**2),
xrange=(-range, range),
yrange=(-range, range),
plot_points=N_sample,fill=True,
contours=[-1.67973392, 1.31869241, 4.31711875,
7.31554508 , 10.31397142, 13.31239775,
16.31082409, 19.30925042],
cmap="jet",
)
plt2=contour_plot(
sign(Bz)*log(Bz**2),
xrange=(-range, range),yrange=(-range, range),
plot_points=N_sample,fill=True,
contours=[
-9.65789873, -6.8984991, -4.13909946,
-1.37969982, 1.37969982, 4.13909946,
6.8984991, 9.65789873],
cmap="jet_r",
)
show(plt2+plt)
v=0.1
gamma=1/sqrt(1-v**2)
t=1
R1=sqrt(gamma**2*(x - v*t)**2+(y + 10^(-6))**2)
R2=sqrt(gamma**2*(x + v*t)**2+(y + 10^(-6))**2)
Ex= (x - v*t)/R1**3 - (x + v*t)/R2**3
Ey= (y+10^(-6))*(1/R1**3 - 1/R2**3)
Bz= v*(y+10^(-6))*(1/R1**3 + 1/R2**3)
for range in (2,20,200,2000):
plt=contour_plot(
log((Ex**2+Ey**2+Bz**2)/2),
xrange=(-range, range),
yrange=(-range, range),
plot_points=N_sample,fill=True,
contours=[-1.67973392, 1.31869241, 4.31711875,
7.31554508 , 10.31397142, 13.31239775,
16.31082409, 19.30925042],
cmap="jet",
)
show(plt)