# 2D-kromme

import numpy as np

E = 0.999
J = 4
Rs = 1

K = E
U = 1 - K**2

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

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)

if n1 < 0:
 n1 = 0.000000001

ubegin = 1.000000001 * n1
ueinde = 0.999999999 * n2
aantalstappen = 1000000
stap = (ueinde - ubegin) / aantalstappen

phi0 = 0

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

for j in range (1, 3, 1):

 for i in range (1, aantalstappen, 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) * np.abs (stap) / 2

  x = r * np.cos (phi)
  y = r * np.sin (phi)

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

 stap = - stap