diff --git a/ffprime/electrostatics/__init__.py b/ffprime/electrostatics/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/ffprime/electrostatics/__pycache__/__init__.cpython-310.pyc b/ffprime/electrostatics/__pycache__/__init__.cpython-310.pyc
new file mode 100644
index 0000000..92e328c
Binary files /dev/null and b/ffprime/electrostatics/__pycache__/__init__.cpython-310.pyc differ
diff --git a/ffprime/electrostatics/__pycache__/multipole.cpython-310.pyc b/ffprime/electrostatics/__pycache__/multipole.cpython-310.pyc
new file mode 100644
index 0000000..3c50ddf
Binary files /dev/null and b/ffprime/electrostatics/__pycache__/multipole.cpython-310.pyc differ
diff --git a/ffprime/electrostatics/multipole.py b/ffprime/electrostatics/multipole.py
new file mode 100644
index 0000000..13311d3
--- /dev/null
+++ b/ffprime/electrostatics/multipole.py
@@ -0,0 +1,20 @@
+import numpy as np
+
+
+def monopole_potential(q: float, r_vec: np.ndarray) -> float:
+    r = np.linalg.norm(r_vec)
+
+    if r < 1e-12:
+        return 0.0
+
+    return q / r
+
+
+def monopole_field(q: float, r_vec: np.ndarray) -> np.ndarray:
+    r = np.linalg.norm(r_vec)
+
+    if r < 1e-12:
+        return np.zeros_like(r_vec)
+
+    return q * r_vec / r**3
+    
\ No newline at end of file
diff --git a/tests/test_electrostatics.py b/tests/test_electrostatics.py
new file mode 100644
index 0000000..7a6130d
--- /dev/null
+++ b/tests/test_electrostatics.py
@@ -0,0 +1,28 @@
+import numpy as np
+from ffprime.electrostatics.multipole import (
+    monopole_potential,
+    monopole_field,
+)
+
+
+def test_monopole_potential_unit_distance():
+    q = 1.0
+    r_vec = np.array([1.0, 0.0, 0.0])
+    V = monopole_potential(q, r_vec)
+    assert np.isclose(V, 1.0)
+
+
+def test_monopole_field_direction():
+    q = 1.0
+    r_vec = np.array([1.0, 0.0, 0.0])
+    E = monopole_field(q, r_vec)
+    assert np.allclose(E, np.array([1.0, 0.0, 0.0]))
+
+
+def test_zero_distance_safe():
+    q = 1.0
+    r_vec = np.array([0.0, 0.0, 0.0])
+    V = monopole_potential(q, r_vec)
+    E = monopole_field(q, r_vec)
+    assert V == 0.0
+    assert np.allclose(E, np.zeros(3))
\ No newline at end of file