# 2D-kromme

import numpy as np

phi0 = 0
E = 0.999
J = 1
K = E
Rs = 1
U = 1 - K**2

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

n1a = 1 / (3 * Rs)
n1b = 1 / (27 * Rs**3) + 1 / (3 * Rs * J**2) - E**2 / (2 * Rs * J**2)
n1c = np.sqrt (E**4 / (4 * Rs**2 * J**4) - E**2 / (3 * Rs**2 * J**4) - E**2 / (27 * Rs**4 * J**2) \
 + 1 / (27 * J**6) + 2 / (27 * Rs**2 * J**4) + 1 / (27 * Rs**4 * J**2))
n1d = n1b + n1c
n1e = n1b - n1c

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

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

n1 = n1a + n1d + n1e

ubegin = 1.000000001 * n1
ueinde = 1.000000001
aantalstappen = 10000000
stap = (ueinde - ubegin) / aantalstappen

u = ubegin
fu = 1 / np.sqrt (a * u**3 + b * u**2 + c * u + d)
phi = phi0

for i in range (1, aantalstappen + 1, 1):

 u = u + stap

 r = 1 / u
 fuprev = fu
 fu = 1 / np.sqrt (a * u**3 + b * u**2 + c * u + d)
 phi = phi + (fuprev + fu) * stap / 2
 x = r * np.cos (phi)
 y = r * np.sin (phi)

 print ('{:10.5e}'.format (x), '{:10.5e}'.format (y))