# 2D-kromme

import numpy as np

def my_arsinh (x):
 y = np.log (x + np.sqrt (x**2 + 1))
 return y

J = 2.1
Rs = 1

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

n1 = (1 - 2 * h2) / (3 * Rs)
n2 = (1 + h2) / (3 * Rs)

if n1 < 0:
 n1 = 0.001

rbegin = 0.999999999 / n1
reinde = 1.000000001 / n2
drempel1 = 1.1 * reinde
drempel2 = 1.01 * reinde
drempel3 = 1.001 * reinde
aantalstappen = 10000000
stap = (reinde - rbegin) / aantalstappen

phi0 = my_arsinh (1)
r = rbegin - stap

while r > reinde:

 if r < drempel1:

  stap = stap / 100
  drempel1 = 0

 if r < drempel2:

  stap = stap / 100
  drempel2 = 0

 if r < drempel3:

  stap = stap / 100
  drempel3 = 0

 r = r + stap
 phi = 2 / h3 * (my_arsinh (h3 / np.sqrt (- Rs / r + 1 / 3 + h2 / 3)) - phi0)
 x = r * np.cos (phi)
 y = r * np.sin (phi)

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