# veldlijnen.py - veld, potentiaal en veldlijnen van twee puntladingen
import numpy as np

K = 8.99e9                                        # 1/(4 pi eps0) in N m²/C²
ladingen = [(+2e-9, np.array([-0.05, 0.0])),     # (lading in C, plaats in m)
            (-2e-9, np.array([+0.05, 0.0]))]


def veld(p):
    """Elektrisch veld (V/m) in punt p: de vectorsom van de bijdragen (superpositie)."""
    E = np.zeros(2)
    for q, a in ladingen:
        r = p - a
        E += K * q * r / np.linalg.norm(r) ** 3
    return E


def potentiaal(p):
    """Potentiaal (V) in punt p, nul ver weg."""
    return sum(K * q / np.linalg.norm(p - a) for q, a in ladingen)


def veldlijn(start, h=0.0005, n_max=3000):
    """Volg de richting van E met Runge-Kutta (vierde orde) in stapjes van h meter."""
    def richting(p):
        E = veld(p)
        return E / np.linalg.norm(E)
    p = np.array(start, float)
    lijn = [p]
    for _ in range(n_max):
        k1 = richting(p)
        k2 = richting(p + h / 2 * k1)
        k3 = richting(p + h / 2 * k2)
        k4 = richting(p + h * k3)
        p = p + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        lijn.append(p)
        if np.linalg.norm(p - ladingen[1][1]) < 0.004:    # aangekomen op de negatieve lading
            return np.array(lijn), True
        if abs(p[0]) > 0.15 or abs(p[1]) > 0.10:          # buiten het beeld
            return np.array(lijn), False
    return np.array(lijn), False


P = np.array([0.0, 0.05])
R = np.array([-0.02, 0.0])
Ex, Ey = veld(P)
print(f"in P: Ex = {Ex:.0f} V/m, Ey = {Ey:.0f} V/m, V = {potentiaal(P):.1f} V")
print(f"in O: Ex = {veld(np.zeros(2))[0]:.0f} V/m")
print(f"in R: Ex = {veld(R)[0]:.0f} V/m, V = {potentiaal(R):.1f} V")
lijnen = []
for hoek in np.radians(np.arange(0, 360, 22.5)):
    start = ladingen[0][1] + 0.005 * np.array([np.cos(hoek), np.sin(hoek)])
    lijnen.append(veldlijn(start))
aan = sum(1 for lijn, op_min in lijnen if op_min)
print(f"{len(lijnen)} veldlijnen: {aan} eindigen in het beeld op -2 nC, {len(lijnen) - aan} verlaten het beeld")