# Hiermee bereken je parameters met hoge nauwkeurigheid

import decimal

from decimal import *
getcontext().prec = 100

J = Decimal ('5')
Rs = Decimal ('1')

T1 = Decimal ('9') * Decimal (Rs) ** Decimal ('2') * Decimal (J) + Decimal (J) ** Decimal ('3')
T2 = (Decimal (J) ** Decimal ('2') - Decimal ('3') * Decimal (Rs) ** Decimal ('2')) ** Decimal ('1.5')
N = Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J)

Ep2 = Decimal ('2') * (Decimal (T1) + Decimal (T2)) / Decimal (N)
Em2 = Decimal ('2') * (Decimal (T1) - Decimal (T2)) / Decimal (N)
Ep = Decimal (Ep2) ** Decimal ('0.5')
Em = Decimal (Em2) ** Decimal ('0.5')

print ("Ep2 = ", Ep2)
print ("Em2 = ", Em2)
print ("Ep = ", Ep)
print ("Em = ", Em)

Dp = (- Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') * Decimal (Ep) ** Decimal ('4') \
 + (Decimal ('36') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 + Decimal ('4') * Decimal (J) ** Decimal ('4')) * Decimal (Ep) ** Decimal ('2') \
 - Decimal ('4') * Decimal (Rs) ** Decimal ('4') \
 - Decimal ('8') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 - Decimal ('4') * Decimal (J) ** Decimal ('4')) / Decimal (J) ** Decimal ('6')

Dm = (- Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') * Decimal (Em) ** Decimal ('4') \
 + (Decimal ('36') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 + Decimal ('4') * Decimal (J) ** Decimal ('4')) * Decimal (Em) ** Decimal ('2') \
 - Decimal ('4') * Decimal (Rs) ** Decimal ('4') \
 - Decimal ('8') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 - Decimal ('4') * Decimal (J) ** Decimal ('4')) / Decimal (J) ** Decimal ('6')

print ("Dp = ", Dp)
print ("Dm = ", Dm)

# Om de discriminant negatief te houden moet Ep groter worden en Em kleiner

Ep = Decimal (Ep) + Decimal ('0.000001')
Em = Decimal (Em) - Decimal ('0.000001')

Dp = (- Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') * Decimal (Ep) ** Decimal ('4') \
 + (Decimal ('36') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 + Decimal ('4') * Decimal (J) ** Decimal ('4')) * Decimal (Ep) ** Decimal ('2') \
 - Decimal ('4') * Decimal (Rs) ** Decimal ('4') \
 - Decimal ('8') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 - Decimal ('4') * Decimal (J) ** Decimal ('4')) / Decimal (J) ** Decimal ('6')

Dm = (- Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') * Decimal (Em) ** Decimal ('4') \
 + (Decimal ('36') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 + Decimal ('4') * Decimal (J) ** Decimal ('4')) * Decimal (Em) ** Decimal ('2') \
 - Decimal ('4') * Decimal (Rs) ** Decimal ('4') \
 - Decimal ('8') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 - Decimal ('4') * Decimal (J) ** Decimal ('4')) / Decimal (J) ** Decimal ('6')

print ("Dp = ", Dp)
print ("Dm = ", Dm)

E = Decimal ('2.01177698018854097539')
J = Decimal ('5')
Rs = Decimal ('1')

D = (- Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') * Decimal (E) ** Decimal ('4') \
 + (Decimal ('36') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 + Decimal ('4') * Decimal (J) ** Decimal ('4')) * Decimal (E) ** Decimal ('2') \
 - Decimal ('4') * Decimal (Rs) ** Decimal ('4') \
 - Decimal ('8') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('2') \
 - Decimal ('4') * Decimal (J) ** Decimal ('4')) / Decimal (J) ** Decimal ('6')

n1a = Decimal ('1') / (Decimal ('3') * Decimal (Rs))
n1b = Decimal ('1') / (Decimal ('27') * Decimal (Rs) ** Decimal ('3')) \
 + Decimal ('1') / (Decimal ('3') * Decimal (Rs) * Decimal (J) ** Decimal ('2')) \
 - Decimal (E) ** Decimal ('2') / (Decimal ('2') * Decimal (Rs) * Decimal (J) ** Decimal ('2'))
n1c = (Decimal (E) ** Decimal ('4') / (Decimal ('4') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('4')) \
 - Decimal (E) ** Decimal ('2') / (Decimal ('3') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('4')) \
 - Decimal (E) ** Decimal ('2') / (Decimal ('27') * Decimal (Rs) ** Decimal ('4') * Decimal (J) ** Decimal ('2')) \
 + Decimal ('1') / (Decimal ('27') * Decimal (J) ** Decimal ('6')) \
 + Decimal ('2') / (Decimal ('27') * Decimal (Rs) ** Decimal ('2') * Decimal (J) ** Decimal ('4')) \
 + Decimal ('1') / (Decimal ('27') * Decimal (Rs) ** Decimal ('4') * Decimal (J) ** Decimal ('2'))) ** Decimal ('0.5')
n1d = Decimal (n1b) + Decimal (n1c)
n1e = Decimal (n1b) - Decimal (n1c)

if n1d < 0:
 n1d = - Decimal (- n1d) ** (Decimal ('1') / Decimal ('3'))
else:
 n1d = Decimal (n1d) ** (Decimal ('1') / Decimal ('3'))

if n1e < 0:
 n1e = - Decimal (- n1e) ** (Decimal ('1') / Decimal ('3'))
else:
 n1e = Decimal (n1e) ** (Decimal ('1') / Decimal ('3'))

n1 = Decimal (n1a) + Decimal (n1d) + Decimal (n1e)

print ("E = ", E)
print ("J = ", J)
print ("Rs = ", Rs)
print ("D = ", D)
print ("n1 = ", n1)