lorentz_transform.py

"""
Lorentz transform of 4-vector field. Reference:
- Introduction to Electrodynamics by Griffiths, chapter 10-12
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import sympy as sp
from sympy.parsing.sympy_parser import parse_expr


def get_transformed_four_potential(*expr, 
                                   vel=sp.Symbol('v'), 
                                   c = sp.Symbol('c')):
    t = sp.Symbol('t')
    x = sp.Symbol('x')
    y = sp.Symbol('y')
    z = sp.Symbol('z')
    if len(expr) < 4:
        expr = expr + tuple(['0' for i in range(4 - len(expr))])
    A = sp.Matrix([[parse_expr(expr[i]) if isinstance(expr[i], str) 
                    else expr[i] for i in range(4)]])
    Λ = sp.Matrix([[1/sp.sqrt(1 - (vel/c)**2), 
                    (vel/c)/sp.sqrt(1 - (vel/c)**2), 0, 0],
                [(vel/c)/sp.sqrt(1 - (vel/c)**2), 
                  1.0/sp.sqrt(1 - (vel/c)**2), 0, 0],
                [0, 0, 1, 0],
                [0, 0, 0, 1]])
    invΛ = Λ**-1
    v = sp.Matrix([[t, x, y, z]])
    v_prime = invΛ*v.T
    t_prime = v_prime[0]
    x_prime = v_prime[1]
    y_prime = v_prime[2]
    z_prime = v_prime[3]
    A_prime = sp.Matrix([[A[i].subs([
                        ('x', x_prime), ('y', y_prime), 
                        ('z', z_prime), ('t', t_prime)
                        ]) for i in range(4)]])
    A_prime = Λ*A_prime.T
    return A_prime[0], A_prime[1], A_prime[2], A_prime[3]


def get_magnetic_field(Ax, Ay, Az):
    dAydz, dAzdy = sp.diff(Ay, 'z'), sp.diff(Az, 'y')
    dAzdx, dAxdz = sp.diff(Az, 'x'), sp.diff(Ax, 'z')
    dAxdy, dAydx = sp.diff(Ax, 'y'), sp.diff(Ay, 'x')
    return -dAydz + dAzdy, -dAzdx + dAxdz, -dAxdy + dAydx


def get_electric_field(V, Ax, Ay, Az):
    dVdx, dVdy, dVdz = sp.diff(V, 'x'), sp.diff(V, 'y'), sp.diff(V, 'z')
    dAxdt, dAydt, dAzdt = sp.diff(Ax, 't'), sp.diff(Ay, 't'), sp.diff(Az, 't')
    return -dVdx - dAxdt, -dVdy - dAydt, -dVdz - dAzdt

def to_func(expr):
    return sp.lambdify(['t', 'x', 'y', 'z'], expr)


potential = 'cos(t - sqrt(x**2 + y**2 + z**2))/sqrt(x**2 + y**2 + z**2)'
potential2 = '(0.5 - 0.5*tanh(-100*(t - sqrt(x**2 + y**2 + z**2))))/sqrt(x**2+y**2+z**2)'
potential3 = '20/sqrt(x**2+y**2+z**2)'
potential4 = '(1.0 + (0.5 - 0.5*tanh(-100*(t - sqrt(x**2 + y**2 + z**2)))))/sqrt(x**2+y**2+z**2)'
V, Ax, Ay, Az = get_transformed_four_potential(potential, vel=0.89, c=1.0)
Ex, Ey, Ez = get_electric_field(V, Ax, Ay, Az)
Bx, By, Bz = get_magnetic_field(Ax, Ay, Az)
V_func, Ax_func, Ay_func, Az_func = [to_func(e) for e in (V, Ax, Ay, Az)]
Ex_func, Ey_func, Ez_func = [to_func(e) for e in (Ex, Ey, Ez)]
Bx_func, By_func, Bz_func = [to_func(e) for e in (Bx, By, Bz)]
x = np.array([np.linspace(-5, 5, 101) for i in range(101)])
y = x.T
x_image = np.array([np.linspace(-5, 5, 101) for i in range(101)])
y_image = x_image.T
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
t = 1.0
z = 1.0
ax.imshow(V_func(t, x_image, y_image, z), extent=(-5.0, 5.0, -5.0, 5.0))
# ax.quiver(x, y, Ax_func(t, x, y, z), Ay_func(t, x, y, z), color='gray', alpha=0.5)
ax.quiver(x, y, Ex_func(t, x, y, z), Ey_func(t, x, y, z), color='red')
# ax.quiver(x, y, Bx_func(t, x, y, z), By_func(t, x, y, z), color='green')
plt.show()
plt.close()


fig = plt.figure()
axes3d = Axes3D(fig)
z = np.array([[np.linspace(-5, 5, 5) 
                for i in range(5)]
                for j in range(5)])
y = np.transpose(z, (0, 2, 1))
x = np.transpose(z, (2, 1, 0))
axes3d.quiver(x, y, z, 
              Ax_func(t, x, y, z), 
              Ay_func(t, x, y, z),
              Az_func(t, x, y, z), color='grey', 
              alpha=0.5)
axes3d.quiver(x, y, z, 
              Ex_func(t, x, y, z), 
              Ey_func(t, x, y, z),
              Ez_func(t, x, y, z), color='red')
axes3d.quiver(x, y, z, 
              Bx_func(t, x, y, z), 
              By_func(t, x, y, z),
              Bz_func(t, x, y, z))
plt.show()
plt.close()







# vel = 0.0
# x_transform = np.array([[1.0/np.sqrt(1.0 - vel**2), vel/np.sqrt(1.0 - vel**2), 0.0, 0.0],
#                         [vel/np.sqrt(1.0 - vel**2), 1.0/np.sqrt(1.0 - vel**2), 0.0, 0.0],
#                         [0.0, 0.0, 1.0, 0.0],
#                         [0.0, 0.0, 0.0, 1.0]])
# invx_transform = np.array([[1.0/np.sqrt(1.0 - vel**2), -vel/np.sqrt(1.0 - vel**2), 0.0, 0.0],
#                           [-vel/np.sqrt(1.0 - vel**2), 1.0/np.sqrt(1.0 - vel**2), 0.0, 0.0],
#                           [0.0, 0.0, 1.0, 0.0],
#                           [0.0, 0.0, 0.0, 1.0]])

# length = 10.0
# n = 11
# z = np.array([[[np.linspace(-length/2, length/2, n) 
#                 for i in range(n)]
#                 for j in range(n)]
#                 for k in range(n)])
# y = np.transpose(z, (0, 1, 3, 2))
# x = np.transpose(y, (0, 2, 1, 3))
# t = np.transpose(x, (1, 0, 2, 3))
# txyz = np.array([t, x, y, z])
# A = np.sin(3*np.pi*t/10.0)
# A = 0.25/np.sqrt(x**2 + y**2 + z**2 + 0.001)
# A = np.array([A, np.zeros([n, n, n, n]), 
#               np.zeros([n, n, n, n]), np.zeros([n, n, n, n])])

# # X_{\mu(ijkl)} T^{\mu}_{\nu} -> X'_{\nu(ijkl)}  
# txyz_prime = np.tensordot(txyz, invx_transform, axes=(0, 0))
# t_prime = txyz_prime[:, :, :, :, 0]
# x_prime = txyz_prime[:, :, :, :, 1]
# y_prime = txyz_prime[:, :, :, :, 2]
# z_prime = txyz_prime[:, :, :, :, 3]
# A_prime = np.sin(3*np.pi*t_prime/10.0)
# # A_prime = 1.0/np.sqrt(x_prime**2 + y_prime**2 + z_prime**2 + 0.5)
# A_prime = np.array([A_prime, np.zeros([n, n, n, n]), 
#                     np.zeros([n, n, n, n]), np.zeros([n, n, n, n])])
# A_prime = np.tensordot(A_prime, x_transform, axes=(0, 0))

# def plot_interact():
#     pass


# fig = plt.figure()
# axes = fig.subplots(1, 2)
# axes[0].set_aspect('equal')
# axes[1].set_aspect('equal')
# for i in range(n):
#     axes[0].plot(txyz[1, i, :, 0, 0], txyz[0, i, :, 0, 0], 
#                 color="black", linestyle='--', linewidth=0.25)
#     axes[0].plot(txyz[1, :, i, 0, 0], txyz[0, :, i, 0, 0],
#                 color="black", linestyle='--', linewidth=0.25)
# axes[0].quiver(x[0, :, 4, 4], t[:, 0, 4, 4], 
#                 A_prime[:, :, 4, 4, 1], A_prime[:, :, 4, 4, 0], scale=4.0)
# axes[0].set_xlabel("x")
# axes[0].set_ylabel("t")
# t_i = 5
# z_i = 5
# axes[1].contourf(x[t_i, :, :, z_i], y[t_i, :, :, z_i], 
#                 A_prime[t_i, :, :, z_i, 0])
# axes[1].quiver(x[t_i, :, :, z_i], y[t_i, :, :, z_i], 
#                 A_prime[t_i, :, :, z_i, 1], A_prime[t_i, :, :, z_i, 2], scale=4.0)
# # axes[1].imshow(A_prime[t_i, :, :, z_i, 0].T, 
# #                 extent=(-5.0, 5.0, -5.0, 5.0))
# axes[1].set_xlabel("x")
# axes[1].set_ylabel("y")
# plt.show()
# plt.close()

# fig = plt.figure()
# axes3d = Axes3D(fig)
# t = 9
# axes3d.quiver(x[0], y[0], z[0], 
#                  A_prime[t, :, :, :, 1], A_prime[t, :, :, :, 2],
#                  A_prime[t, :, :, :, 3], length=4
#                    )
# plt.show()
# plt.close()

lorentz_transform2d.py

import numpy as np
import matplotlib.pyplot as plt


vel = 0.5
lorentz_transform = np.array([[1.0/np.sqrt(1.0 - vel**2), vel/np.sqrt(1.0 - vel**2)],
                              [vel/np.sqrt(1.0 - vel**2), 1.0/np.sqrt(1.0 - vel**2)]])
inv_lorentz_transform = np.array([[1.0/np.sqrt(1.0 - vel**2), -vel/np.sqrt(1.0 - vel**2)],
                                  [-vel/np.sqrt(1.0 - vel**2), 1.0/np.sqrt(1.0 - vel**2)]])

n = 11
x = np.array([np.linspace(-5.0, 5.0, n) for _ in range(n)])
y = x.T.copy()
r = np.array([x, y])
# (i, j, k) (i, l)
r_prime = np.tensordot(r, inv_lorentz_transform, axes=(0, 1))
x_prime = r_prime[:,:,0]
y_prime = r_prime[:,:,1]
# for line in zip(x_prime, y_prime):
#     plt.plot(*line)
u = np.zeros([n, n])
v = np.array([[1.0/abs(x_i) if x_i != 0 else 2.0 for x_i in row] for row in x_prime])
uv = np.array([u, v])
uv_prime = np.tensordot(uv, lorentz_transform, axes=(0, 1))
u_prime = uv_prime[:,:,0]
v_prime = uv_prime[:,:,1]
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.quiver(x, y, u_prime, v_prime)
plt.grid()
plt.show()
plt.close()