https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates

import numpy as np
import matplotlib.pyplot as plt

# %matplotlib widget

Field around a point charge

q1 = np.array([0, 0, 0])
A = np.array([-7, -5, 6])

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))
#point charge
ax.plot(*q1, color = 'red', marker = 'o')
ax.text(q1[0] + 0.25, q1[1] + 0.25, q1[2] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

radius_vector = A - q1
radius_scalar = np.linalg.norm(radius_vector)

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))
#point charge
ax.plot(*q1, color = 'red', marker = 'o')
ax.text(q1[0] + 0.25, q1[1] + 0.25, q1[2] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

#radius vector
ax.quiver(*q1, *radius_vector, color = 'black')

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

field_unit_vector = radius_vector / radius_scalar

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))
#point charge
ax.plot(*q1, color = 'red', marker = 'o')
ax.text(q1[0] + 0.25, q1[1] + 0.25, q1[2] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

#radius vector
ax.quiver(*q1, *radius_vector, color = 'black')

#field vector
ax.quiver(*A, *field_unit_vector, color = 'red')

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

theta = np.linspace(0, 2*np.pi, 100)
alpha = np.linspace(0, np.pi, 100)
X = q1[0] + radius_scalar * np.outer(np.cos(theta), np.sin(alpha))
Y = q1[1] + radius_scalar * np.outer(np.sin(theta), np.sin(alpha))
Z = q1[2] + radius_scalar * np.outer(np.ones(theta.shape), np.cos(alpha))

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))
#point charge
ax.plot(*q1, color = 'red', marker = 'o')
ax.text(q1[0] + 0.25, q1[1] + 0.25, q1[2] + 0.25, s = '$q_1$', fontsize = 12)
#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)
#radius vector
ax.quiver(*q1, *radius_vector, color = 'black')
#field vector
ax.quiver(*A, *field_unit_vector, color = 'red')
#field sphere
ax.plot_surface(X, Y, Z, color = 'blue', alpha = 0.2)
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

k = 1 / (4*np.pi*8.8541878128e-12)
q = 1.60217662e-19
electric_field = k * q / radius_scalar**2
print(electric_field)

1.3090586684181768e-11

Field around a infinitly line charge

length = 20
center_of_charge = np.array([0, 0, 0])
line_q1 = np.array([[center_of_charge[0], center_of_charge[0], center_of_charge[0]],
                    [center_of_charge[1], center_of_charge[1], center_of_charge[1]],
                    [center_of_charge[2] - (length/2), center_of_charge[2], center_of_charge[2] + (length/2)]])
A = np.array([6, -5, 4])

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))

#line charge
ax.plot(*line_q1, color = 'red')
ax.text(line_q1[0, 1] + 0.25, line_q1[1, 1] + 0.25, line_q1[2, 1] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

A_prime = np.array([center_of_charge[0], center_of_charge[1], center_of_charge[2] + (A[2] -  center_of_charge[2])])
radius_vector = A - A_prime
radius_scalar = np.linalg.norm(radius_vector)

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))

#line charge
ax.plot(*line_q1, color = 'red')
ax.text(line_q1[0, 1] + 0.25, line_q1[1, 1] + 0.25, line_q1[2, 1] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

#radius vector
ax.quiver(*A_prime, *radius_vector, color = 'black')

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

field_unit_vector = radius_vector / radius_scalar

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))

#line charge
ax.plot(*line_q1, color = 'red')
ax.text(line_q1[0, 1] + 0.25, line_q1[1, 1] + 0.25, line_q1[2, 1] + 0.25, s = '$q_1$', fontsize = 12)

#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)

#radius vector
ax.quiver(*A_prime, *radius_vector, color = 'black')

#field vector
ax.quiver(*A, *field_unit_vector, color = 'red')

ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

z = np.linspace(-length/2, length/2, 50)
angle = np.linspace(0, 2*np.pi, 50)
theta, Z = np.meshgrid(angle, z)

X = radius_scalar*np.cos(theta) + center_of_charge[0]
Y = radius_scalar*np.sin(theta) + center_of_charge[1]
Z = Z + center_of_charge[2]

fig = plt.figure(figsize = (6, 6))
ax = plt.axes(projection='3d', box_aspect = (1, 1, 1))
#line charge
ax.plot(*line_q1, color = 'red')
ax.text(line_q1[0, 1] + 0.25, line_q1[1, 1] + 0.25, line_q1[2, 1] + 0.25, s = '$q_1$', fontsize = 12)
#point A
ax.plot(*A, color = 'green', marker = 'o')
ax.text(A[0] + 0.25, A[1] + 0.25, A[2] + 0.25, s = '$A$', fontsize = 12)
#radius vector
ax.quiver(*A_prime, *radius_vector, color = 'black')
#field vector
ax.quiver(*A, *field_unit_vector, color = 'red')
#field cylinder
ax.plot_surface(X, Y, Z, color = 'blue', alpha = 0.2)
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_zlim(-10, 10)
plt.grid()

k = 1 / (2*np.pi*8.8541878128e-12)
q = 1.60217662e-19
electric_field = k * q / (radius_scalar*length) 

print(electric_field)

1.1999704460499954e-10