# mpp.py - het punt van het grootste vermogen van een zonnepaneel, drie numerieke methoden
import math

ISC, VOC, A = 10.9, 41.0, 2.44                 # kortsluitstroom (A), openklemspanning (V), diodeparameter (V)
I0 = ISC / (math.exp(VOC / A) - 1)             # verzadigingsstroom: I(VOC) = 0


def P(u):                                      # vermogen (W) bij spanning u (V)
    return u * (ISC - I0 * (math.exp(u / A) - 1))


def dP(u, h=1e-5):                             # afgeleide: centrale differentie (hoofdstuk 1)
    return (P(u + h) - P(u - h)) / (2 * h)


# 1. bisectie op de afgeleide: dP > 0 links van het MPP, dP < 0 rechts
a, b, n = 0.0, VOC, 0
while b - a > 0.001:
    m = (a + b) / 2
    if dP(m) > 0:
        a = m
    else:
        b = m
    n += 1
print("bisectie op P':  U = %.3f V, P = %.2f W na %d stappen" % ((a + b) / 2, P((a + b) / 2), n))

# 2. gulden snede op P zelf (geen afgeleide nodig)
g = (math.sqrt(5) - 1) / 2
a, b, n = 0.0, VOC, 0
x1, x2 = b - g * (b - a), a + g * (b - a)
p1, p2 = P(x1), P(x2)
while b - a > 0.01:
    if p1 > p2:
        b, x2, p2 = x2, x1, p1
        x1 = b - g * (b - a)
        p1 = P(x1)
    else:
        a, x1, p1 = x1, x2, p2
        x2 = a + g * (b - a)
        p2 = P(x2)
    n += 1
print("gulden snede:     U = %.3f V, P = %.2f W na %d stappen" % ((a + b) / 2, P((a + b) / 2), n))

# 3. perturb and observe, zoals een MPPT-regelaar: kleine stap, kijk of P stijgt
u, stap, richting = 30.0, 0.5, +1
vorig = P(u)
spoor = []
for k in range(14):
    u += richting * stap
    nu = P(u)
    if nu < vorig:                             # slechter: keer om
        richting = -richting
    vorig = nu
    spoor.append("%.1f" % u)
print("perturb and observe:", " ".join(spoor))