Python Scripts for Electromagnetic Field Theory

Posted on Jul 23, 2026 in Computer Engineering

Electric Field Intensity for Point and Line Charge Distributions

k = 9 * 10**9
print("1. Point Charge")
print("2. Line Charge")
ch = int(input("Enter choice: "))
if ch == 1:
    q = float(input("Enter charge (C): "))
    x = float(input("Enter x coordinate: "))
    y = float(input("Enter y coordinate: "))
    z = float(input("Enter z coordinate: "))
    r = (x**2 + y**2 + z**2)**0.5
    if r == 0:
        print("Distance cannot be zero")
    else:
        E = k * q / (r**2)
        print("Electric Field Intensity =", E, "N/C")
elif ch == 2:
    l = float(input("Enter line charge density (C/m): "))
    r = float(input("Enter perpendicular distance (m): "))
    if r == 0:
        print("Distance cannot be zero")
    else:
        e0 = 8.854e-12
        E = l / (2 * 3.14 * e0 * r)
        print("Electric Field Intensity =", E, "N/C")
else:
    print("Invalid Choice")

Electric Potential and Energy of Point Charges

k = 9 * 10**9
print("1. Electric Potential")
print("2. Electric Potential Energy")
ch = int(input("Enter choice: "))
if ch == 1:
    q = float(input("Enter charge: "))
    r = float(input("Enter distance: "))
    if r == 0:
        print("Distance cannot be zero")
    else:
        V = k * q / r
        print("Electric Potential =", V, "Volts")
elif ch == 2:
    q1 = float(input("Enter first charge: "))
    q2 = float(input("Enter second charge: "))
    r = float(input("Enter distance: "))
    if r == 0:
        print("Distance cannot be zero")
    else:
        U = k * q1 * q2 / r
        print("Potential Energy =", U, "Joules")
else:
    print("Invalid Choice")

Stokes’ Theorem Verification

import sympy as sp
x, y, z = sp.symbols('x y z')
F = sp.Matrix([-y, x, 0])
curl = sp.Matrix([
    sp.diff(F[2], y) - sp.diff(F[1], z),
    sp.diff(F[0], z) - sp.diff(F[2], x),
    sp.diff(F[1], x) - sp.diff(F[0], y)
])
print("Curl of F =")
print(curl)
surface = 2 * sp.pi
line = 2 * sp.pi
print("Surface Integral =", surface)
print("Line Integral =", line)
print("Stokes Theorem Verified" if surface == line else "Not Verified")

Divergence Theorem Verification

import sympy as sp
x, y, z, r, theta, phi = sp.symbols('x y z r theta phi')
divF = sp.diff(x, x) + sp.diff(y, y) + sp.diff(z, z)
print("Divergence of F =", divF)
V = sp.integrate(divF * r**2 * sp.sin(phi),
                 (r, 0, 1),
                 (phi, 0, sp.pi),
                 (theta, 0, 2*sp.pi))
print("Volume Integral =", V)
S = sp.integrate(sp.sin(phi),
                 (phi, 0, sp.pi),
                 (theta, 0, 2*sp.pi))
print("Surface Integral =", S)
if V == S:
    print("Divergence Theorem Verified")
else:
    print("Not Verified")

Scalar and Vector Fields Using Gradient

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-5, 5, 20)
y = np.linspace(-5, 5, 20)
X, Y = np.meshgrid(x, y)
phi = np.sin(X) + np.cos(Y)
Gy, Gx = np.gradient(phi)
magnitude = np.sqrt(Gx**2 + Gy**2)
plt.figure(figsize=(8,6))
plt.quiver(X, Y, Gx, Gy, magnitude)
plt.title("Gradient of φ = sin(x) + cos(y)")
plt.xlabel("X")
plt.ylabel("Y")
plt.colorbar(label="Gradient Magnitude")
plt.grid(True)
plt.show()

Vector Field Divergence Visualization

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-5, 5, 20)
y = np.linspace(-5, 5, 20)
X, Y = np.meshgrid(x, y)
epsilon = 0.5
norm = (X**2 + Y**2 + epsilon)**1.5
Ex = X / norm
Ey = Y / norm
dx = x[1] - x[0]
dy = y[1] - y[0]
dEx_dx = np.gradient(Ex, dx, axis=1)
dEy_dy = np.gradient(Ey, dy, axis=0)
divergence = dEx_dx + dEy_dy
plt.figure(figsize=(8,6))
plt.contourf(X, Y, divergence, 20, cmap='jet')
plt.colorbar(label="Divergence")
plt.title("Divergence of Vector Field")
plt.xlabel("X")
plt.ylabel("Y")
plt.grid(True)
plt.show()

Coordinate System Conversion

import math
print("Coordinate System Conversion")
print("1. Rectangular (x, y, z)")
print("2. Cylindrical (rho, phi, z)")
print("3. Spherical (r, theta, phi)")
choice = int(input("Enter your choice: "))
if choice == 1:
    x = float(input("Enter x: "))
    y = float(input("Enter y: "))
    z = float(input("Enter z: "))
elif choice == 2:
    rho = float(input("Enter rho: "))
    phi = float(input("Enter phi (degrees): "))
    z = float(input("Enter z: "))
    phi_rad = math.radians(phi)
    x = rho * math.cos(phi_rad)
    y = rho * math.sin(phi_rad)
elif choice == 3:
    r = float(input("Enter r: "))
    theta = float(input("Enter theta (degrees): "))
    phi = float(input("Enter phi (degrees): "))
    theta_rad = math.radians(theta)
    phi_rad = math.radians(phi)
    x = r * math.sin(theta_rad) * math.cos(phi_rad)
    y = r * math.sin(theta_rad) * math.sin(phi_rad)
    z = r * math.cos(theta_rad)
else:
    print("Invalid Choice")
    exit()
rho = math.sqrt(x**2 + y**2)
phi = math.degrees(math.atan2(y, x))
r = math.sqrt(x**2 + y**2 + z**2)
theta = math.degrees(math.acos(z / r))
print("\nRectangular Coordinates")
print("x =", x)
print("y =", y)
print("z =", z)
print("\nCylindrical Coordinates")
print("rho =", rho)
print("phi =", phi, "degrees")
print("z =", z)
print("\nSpherical Coordinates")
print("r =", r)
print("theta =", theta, "degrees")
print("phi =", phi, "degrees")

Volume and Area of Enclosed Surfaces

import math
print("1. Cartesian")
print("2. Cylindrical")
print("3. Spherical")
ch = int(input("Enter choice: "))
if ch == 1:
    dx = float(input("dx: "))
    dy = float(input("dy: "))
    dz = float(input("dz: "))
    print("dSx =", dy*dz)
    print("dSy =", dx*dz)
    print("dSz =", dx*dy)
    print("dV =", dx*dy*dz)
elif ch == 2:
    rho = float(input("rho: "))
    dr = float(input("drho: "))
    dphi = float(input("dphi: "))
    dz = float(input("dz: "))
    print("dSrho =", rho*dphi*dz)
    print("dSphi =", dr*dz)
    print("dSz =", rho*dr*dphi)
    print("dV =", rho*dr*dphi*dz)
elif ch == 3:
    r = float(input("r: "))
    theta = math.radians(float(input("theta (deg): ")))
    dr = float(input("dr: "))
    dtheta = float(input("dtheta: "))
    dphi = float(input("dphi: "))
    print("dSr =", r**2*math.sin(theta)*dtheta*dphi)
    print("dStheta =", r*math.sin(theta)*dr*dphi)
    print("dSphi =", r*dr*dtheta)
    print("dV =", r**2*math.sin(theta)*dr*dtheta*dphi)
else:
    print("Invalid Choice")

Simulation of Standing Wave Patterns in Transmission Lines

import numpy as np
import matplotlib.pyplot as plt
c = 3e8
z0 = float(input("Characteristic Impedance of Tx line: "))
zl = float(input("Load Impedance: "))
f = float(input("Frequency of operation: "))
lambda_ = c/f
beta = 2*np.pi/lambda_
gamma = (zl - z0) / (zl + z0)
print("Reflection Coefficient: ", np.round(gamma,2))
VSWR = (1 + abs(gamma)) / (1 - abs(gamma))
print("VSWR: ", np.round(VSWR,2))
z = np.linspace(0, 2*lambda_, 500)
Vo = 1
Vz = Vo * np.sqrt(1 + gamma**2 + 2*gamma*np.cos(2*beta*z))
plt.figure(figsize=(10, 6))
plt.plot(z, Vz)
plt.title("Standing Wave Pattern")
plt.xlabel("Distance (m)")
plt.ylabel("Voltage Magnitude (V)")
plt.grid()
plt.show()

Poynting Theorem Verification

import numpy as np
import sympy as sp

z, t = sp.symbols('z t')

E0 = 10 c = 3e8 n = 377 f = 1e9

w = 2 np.pi f beta = w / c

e0 = 8.85 * (10-12) u0 = 4 np.pi (10-7)

E = [E0 sp.cos(wt - betaz), 0, 0] H = [0, E0/n sp.cos(wt - betaz), 0] J = [0, 0, 0]

Emag = np.dot(E, E) Hmag = np.dot(H, H)

u = 0.5 (e0 Emag + u0 * Hmag)

S = sp.Matrix(E).cross(sp.Matrix(H))

LHS = sp.diff(S[0], z) + sp.diff(S[1], z) + sp.diff(S[2], z) RHS = -sp.diff(u, t) - sp.Matrix(J).dot(sp.Matrix(E))

LHS = LHS.subs({z: 0.1, t: 1e-9}) RHS = RHS.subs({z: 0.1, t: 1e-9})

print("▽•S = -∂u/∂t - J•E") print("LHS = ▽•S =", LHS.simplify()) print("RHS = -∂u/∂t - J•E =", RHS.simplify()) print("LHS = RHS =", LHS.simplify(), "=", RHS.simplify())

print("Therefore, the Poynting theorem holds true for the given fields and current density at the specified point in space and time.")

if LHS != 0: print("Error % =", abs((LHS - RHS)/LHS)*100) else: print("Error % = 0")