sim.py

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import interp1d

import ipdb

#ref The Finite-Difference Time Domain Method for Solving Maxwell’s Equations

size = 15
dx = 0.1
dy = 0.1
dz = 0.1
d = [dx, dy, dz]
dt = 0.8/(3e8*np.sqrt(1/dx**2+1/dy**2+1/dz**2))
ctr = int(size/2)

X, Y, Z = np.meshgrid(range(size-1), range(size-1), range(size-1), indexing='ij')

# at vertices  x       y     z
Ex = np.zeros([size-1, size, size])
Ey = np.zeros([size, size-1, size])
Ez = np.zeros([size, size, size-1])
E = [Ex, Ey, Ez]
Ep = [f[:] for f in E]
Epp = [f[:] for f in E]

# at faces     x       y       z
Hx = np.zeros([size, size-1, size-1])
Hy = np.zeros([size-1, size, size-1])
Hz = np.zeros([size-1, size-1, size])
H = [Hx, Hy, Hz]

permitivity = np.zeros([size, size, size])+8.85e-12 # epsilon
conductivity = np.zeros([size, size, size]) # sigma
permeability = np.zeros([size, size, size])+np.pi*4e-7 # mu

# coper wire
#conductivity[:,5,5] = 5.96e7
#permeability[:,5,5] = 1.256629e-6

# voltage source
E[0][:,ctr,ctr] = 2

def avg(a, axis):
    "compute midpoints along a given axis"
    idx = [slice(None) for _ in range(a.ndim)]
    idx[axis] = slice(1, None)
    op1 = a[tuple(idx)]
    idx[axis] = slice(None, -1)
    op2 = a[tuple(idx)]
    return (op1+op2)/2

def curl(Fx, Fy, Fz):
    "3D cartesian curl"
    Fzdy = np.diff(Fz, axis=1)/dy
    Fydz = np.diff(Fy, axis=2)/dz
    Fxdz = np.diff(Fx, axis=2)/dz
    Fzdx = np.diff(Fz, axis=0)/dx
    Fydx = np.diff(Fy, axis=0)/dx
    Fxdy = np.diff(Fx, axis=1)/dy
    return [Fzdy-Fydz,
            Fxdz-Fzdx,
            Fydx-Fxdy]

def boundary(E, Ep, Epp):
    "Update boundaries in E using Liao's absorbing boundary condition"
    for i in range(3):
        coord = np.arange(E.shape[i])/d[i]
        ip = interp1d(coord, Ep, kind='quadratic', axis=i)
        ipp = interp1d(coord, Epp, kind='quadratic', axis=i)
        Ec = ip(3e8*dt)
        E2c = ipp(2*3e8*dt)
        idx = [slice(None), slice(None), slice(None)]
        idx[i] = 0
        E[tuple(idx)] = 2*Ec-E2c
        Ec = ip(coord[-1]-3e8*dt)
        E2c = ipp(coord[-1]-2*3e8*dt)
        idx = [slice(None), slice(None), slice(None)]
        idx[i] = -1
        E[tuple(idx)] = 2*Ec-E2c


for t in range(1000):
    Epp = Ep
    Ep = [f[:] for f in E]
    CurlH = curl(H[0][1:-1,:,:], H[1][:,1:-1,:], H[2][:,:,1:-1])
    for i in range(3):
        idx = [slice(1, -1), slice(1, -1), slice(1,-1)]
        idx[i] = slice(None)
        idx = tuple(idx)
        pmt = avg(permitivity, i)[idx]
        cnd = avg(conductivity, i)[idx]
        coef1 = (1-cnd*dt/(2*pmt))/(1+cnd*dt/(2*pmt))
        coef2 = (dt/pmt)/(1+cnd*dt/(2*pmt))
        E[i][idx] = coef1*E[i][idx]+coef2*CurlH[i]

        boundary(E[i], Ep[i], Epp[i])


    E[0][:,ctr,ctr] = 2

    CurlE = curl(E[0], E[1], E[2])
    for i in range(3):
        mu = avg(avg(permeability, (i+1)%3),(i+2)%3)
        H[i] = H[i]-dt/mu*CurlE[i]

    print(*[np.count_nonzero(m) for m in E + H])

if True:
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.quiver(X, Y, Z, avg(H[0], 0), avg(H[1], 1), avg(H[2], 2), normalize=False)

if True:
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.quiver(X, Y, Z, avg(avg(E[0], 1), 2), avg(avg(E[1], 0), 2), avg(avg(E[2], 1), 0), normalize=False)

plt.show()