import numpy as np
import matplotlib.pyplot as plt

P0 = np.array([1, 1, 1])
n = np.array([1, 0, 1])

# d = -(ax_0 + by_0 + cz_0)
d = -np.dot(n, P0)

num_points = 3000
x = np.random.uniform(low=-3, high=3, size=num_points)
y = np.random.uniform(low=-3, high=3, size=num_points)

# ax + by + cz + d = 0
# cz = -(ax + by + d)
# z = -(ax + by + d) / c
a, b, c = n
z = -(a * x + b * y + d) / c

P = np.hstack([x[None, :].T, y[None, :].T, z[None, :].T])
p = P - P0

assert np.allclose(np.dot(p, n), 0)

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection="3d")
ax.scatter(
    *P0, 
    c="red", 
    s=20, 
    marker="o", 
    label=f"P0={P0}"
)

ax.scatter(
    P[:, 0], 
    P[:, 1], 
    P[:, 2], 
    c="blue", 
    s=3, 
    alpha=0.25,
    marker="o", 
    label="P"
)

ax.quiver(
    P0[0], 
    P0[1], 
    P0[2], 
    n[0], 
    n[1], 
    n[2], 
    color="red", 
    arrow_length_ratio=0.2, 
    linewidth=2, 
    label=f"n={n}"
)

ax.set_proj_type("ortho")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title(f"\n Plane Equation: {n[0]}x + {n[1]}y + {n[2]}z + {d} = 0 \n")
ax.legend()
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_zlim(-4, 4)
ax.axis("off")

plt.tight_layout()
plt.show()