# Hiermee bereken je de parameters in het geval dat D > 0

import numpy as np

E = 0.994882114
J = 5
Rs = 1

K = E
U = 1 - K**2

a = Rs
b = -1
c = Rs / J**2
d = - U / J**2

D = 18 * a * b * c * d - 4 * a * c**3 - 27 * a**2 * d**2 + b**2 * c**2 - 4 * b**3 * d

h1 = 1 - 3 * Rs**2 / J**2
h2 = np.sqrt (h1)

T = - 2 * J**2 - 18 * Rs**2 + 27 * Rs**2 * E**2
N = 2 * (J**2 - 3 * Rs**2) * h2

n1 = 1 / (3 * Rs) - 2 / (3 * Rs) * h2 * np.sin ((np.arcsin (T / N) + 1 * np.pi) / 3)
n2 = 1 / (3 * Rs) - 2 / (3 * Rs) * h2 * np.sin ((np.arcsin (T / N) + 3 * np.pi) / 3)
n3 = 1 / (3 * Rs) - 2 / (3 * Rs) * h2 * np.sin ((np.arcsin (T / N) + 5 * np.pi) / 3)

r1 = 1 / n1
r2 = 1 / n2
r3 = 1 / n3

print ("E = ", E)
print ("J = ", J)
print ("Rs = ", Rs)
print ("K = ", K)
print ("U = ", U)
print ("a = ", a)
print ("b = ", b)
print ("c = ", c)
print ("d = ", d)
print ("D = ", D)
print ("n1 = ", n1)
print ("n2 = ", n2)
print ("n3 = ", n3)
print ("r1 = ", r1)
print ("r2 = ", r2)
print ("r3 = ", r3)