電磁場のベクトル表現:直線運動を行う点電荷¶

３次元空間の設定¶

In [1]:
%display latex #latexを使って、結果を表示します。
load( "setupSphericalCoordinates.sage")
\(\displaystyle \verb|SageMath|\verb| |\verb|version|\verb| |\verb|10.0,|\verb| |\verb|Release|\verb| |\verb|Date:|\verb| |\verb|2023-05-20|\)
In [2]:
E.<x,y,z> = EuclideanSpace(name="E")
show(E.atlas(),E.frames(),[type(f) for f in E.frames()])
\(\displaystyle \left[\left(E,(x, y, z)\right)\right] \left[\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right)\right] \left[\verb|<class|\verb| |\verb|'sage.manifolds.differentiable.vectorframe.CoordFrame_with_category'>|\right]\)
In [3]:
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
In [4]:
E.default_frame(),E.default_chart()
Out[4]:
\(\displaystyle \left(\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right), \left(E,(x, y, z)\right)\right)\)
In [5]:
g[NF,:], g.up(g).up(g)[NF,:]
Out[5]:
\(\displaystyle \left(\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right), \left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)\right)\)
In [6]:
CC.frame(),E.frames()
Out[6]:
\(\displaystyle \left(\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right), \left[\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right), \left(E, \left(\frac{\partial}{\partial r },\frac{\partial}{\partial {\theta} },\frac{\partial}{\partial {\phi} }\right)\right), \left(E, \left(e_{ r },e_{ {\theta} },e_{ {\phi} }\right)\right)\right]\right)\)
In [7]:
E.changes_of_frame()
Out[7]:
\(\displaystyle \left\{\left(\left(E, \left(\frac{\partial}{\partial r },\frac{\partial}{\partial {\theta} },\frac{\partial}{\partial {\phi} }\right)\right), \left(E, \left(e_{ r },e_{ {\theta} },e_{ {\phi} }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}, \left(\left(E, \left(e_{ r },e_{ {\theta} },e_{ {\phi} }\right)\right), \left(E, \left(\frac{\partial}{\partial r },\frac{\partial}{\partial {\theta} },\frac{\partial}{\partial {\phi} }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}, \left(\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right), \left(E, \left(\frac{\partial}{\partial r },\frac{\partial}{\partial {\theta} },\frac{\partial}{\partial {\phi} }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}, \left(\left(E, \left(\frac{\partial}{\partial r },\frac{\partial}{\partial {\theta} },\frac{\partial}{\partial {\phi} }\right)\right), \left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}, \left(\left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right), \left(E, \left(e_{ r },e_{ {\theta} },e_{ {\phi} }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}, \left(\left(E, \left(e_{ r },e_{ {\theta} },e_{ {\phi} }\right)\right), \left(E, \left(e_{ x },e_{ y },e_{ z }\right)\right)\right) : \mbox{Field of tangent-space automorphisms on the Euclidean space E}\right\}\)
In [8]:
E.change_of_frame(SPF,CF)[:],E.change_of_frame(CF,SPF)[:],
Out[8]:
\(\displaystyle \left(\left(\begin{array}{rrr} \frac{x}{\sqrt{x^{2} + y^{2} + z^{2}}} & \frac{y}{\sqrt{x^{2} + y^{2} + z^{2}}} & \frac{z}{\sqrt{x^{2} + y^{2} + z^{2}}} \\ \frac{x z}{\sqrt{x^{2} + y^{2} + z^{2}} \sqrt{x^{2} + y^{2}}} & \frac{y z}{\sqrt{x^{2} + y^{2} + z^{2}} \sqrt{x^{2} + y^{2}}} & -\frac{\sqrt{x^{2} + y^{2}}}{\sqrt{x^{2} + y^{2} + z^{2}}} \\ -\frac{y}{\sqrt{x^{2} + y^{2}}} & \frac{x}{\sqrt{x^{2} + y^{2}}} & 0 \end{array}\right), \left(\begin{array}{rrr} \frac{x}{\sqrt{x^{2} + y^{2} + z^{2}}} & \frac{x z}{\sqrt{x^{2} + y^{2} + z^{2}} \sqrt{x^{2} + y^{2}}} & -\frac{y}{\sqrt{x^{2} + y^{2}}} \\ \frac{y}{\sqrt{x^{2} + y^{2} + z^{2}}} & \frac{y z}{\sqrt{x^{2} + y^{2} + z^{2}} \sqrt{x^{2} + y^{2}}} & \frac{x}{\sqrt{x^{2} + y^{2}}} \\ \frac{z}{\sqrt{x^{2} + y^{2} + z^{2}}} & -\frac{\sqrt{x^{2} + y^{2}}}{\sqrt{x^{2} + y^{2} + z^{2}}} & 0 \end{array}\right)\right)\)

等速度運動する点電荷の作る電場¶

In [9]:
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()
In [10]:
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)
In [11]:
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)
In [12]:
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$ では下の図のように、電場強度の分布は進行方向に対して、扁平に圧縮された分布と なっていることがわかります。

In [13]:
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)
In [14]:
#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)
In [15]:
#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)
In [16]:
#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)
In [17]:
#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)
In [18]:
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)
In [19]:
anim=animate(plts)
anim.show(delay=50, iterations=5)
In [20]:
anim.save("./movie/movingCharges.mp4",show_path=True)

Animation saved to file /Users/noboru/src/JupyterNotebooks/SageMath/SageManifolds/movie/movingCharges.mp4.
In [21]:
anim.show(delay=10, iterations=5)
In [24]:
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)
In [25]:
anim=animate(plts)
anim.save("./movie/SlowMovingCharges.mp4",show_path=True)

Animation saved to file /Users/noboru/src/JupyterNotebooks/SageMath/SageManifolds/movie/SlowMovingCharges.mp4.
In [26]:
anim.show(delay=10, iterations=5)
In [27]:
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)
In [28]:
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)