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

# Define lattice constants
a = 1  # lattice constant

# Define Miller indices for the planes
hkl1 = np.array([1, 2, 1])  # (121)
hkl2 = np.array([1, -2, -1])  # (1-2-1)

# Define plane normal vectors
n1 = np.array([hkl1[0], hkl1[1], hkl1[2], 0])
n2 = np.array([hkl2[0], hkl2[1], hkl2[2], 0])

# Define a point on each plane
p1 = np.array([0, 0, 0])
p2 = np.array([a/2, a/2, 0])

# Define function to plot planes
def plot_plane(n, p, color):
    d = -np.dot(n, p)
    xx, yy = np.meshgrid(range(-10, 10), range(-10, 10))
    zz = (-n[0]*xx - n[1]*yy - d)*1. /n[2]
    ax.plot_surface(xx, yy, zz, color=color, alpha=0.3)

# Create 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot cubic lattice
for i in range(-10, 10):
    for j in range(-10, 10):
        for k in range(-10, 10):
            x = i * a
            y = j * a
            z = k * a
            ax.scatter(x, y, z, color='black', marker='o', s=10)

# Plot planes
plot_plane(n1, p1, 'blue')
plot_plane(n2, p2, 'red')

# Set axis labels and limits
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)

# Show plot
plt.show()