
=== ORBIT I (vertex (a,0)), every ellipse ===
  check all_on_ellipse: True
  check closes_after_N: True
  check direction_returns_parallel: True
  check no_earlier_return: True
  check all_sides_tangent: True
  contact ratios: [b/(a + b), 1/2, a/(a + b), b/(a + b), 1/2, a/(a + b)]
  vertices:
    P1 = (a, 0)
    P2 = (a**2/(a + b), b**(3/2)*sqrt(2*a + b)/(a + b))
    P3 = (-a**2/(a + b), b**(3/2)*sqrt(2*a + b)/(a + b))
    P4 = (-a, 0)
    P5 = (-a**2/(a + b), -b**(3/2)*sqrt(2*a + b)/(a + b))
    P6 = (a**2/(a + b), -b**(3/2)*sqrt(2*a + b)/(a + b))
  A = 2*a*b**(3/2)*(2*a + b)**(3/2)/(a + b)**2
  A' = 2*a*sqrt(b)*(a + 2*b)/sqrt(2*a + b)
  A'' = 2*a**2*b**(3/2)*(a + 2*b)**2*sqrt(2*a + b)/(a + b)**4
  A*_O = 2*a*(a**3 - 2*a**2*b + 2*a*b**2 + 2*b**3)/(b**(3/2)*sqrt(2*a + b))
  A'*_O = 2*sqrt(b)*(a + b)**2*(4*a**2 - 2*a*b + b**2)/(a*(2*a + b)**(3/2))
  Aped_O = 2*a*b**(3/2)*(2*a + b)**(3/2)*(a**2 + a*b + b**2)/(a + b)**4
  A'ped_O = sqrt(b)*(2*a + b)**(3/2)/2

=== ORBIT II (vertex (0,b)), every ellipse ===
  check all_on_ellipse: True
  check closes_after_N: True
  check direction_returns_parallel: True
  check no_earlier_return: True
  check all_sides_tangent: True
  contact ratios: [a/(a + b), 1/2, b/(a + b), a/(a + b), 1/2, b/(a + b)]
  vertices:
    P1 = (0, b)
    P2 = (-a**(3/2)*sqrt(a + 2*b)/(a + b), b**2/(a + b))
    P3 = (-a**(3/2)*sqrt(a + 2*b)/(a + b), -b**2/(a + b))
    P4 = (0, -b)
    P5 = (a**(3/2)*sqrt(a + 2*b)/(a + b), -b**2/(a + b))
    P6 = (a**(3/2)*sqrt(a + 2*b)/(a + b), b**2/(a + b))
  A = 2*a**(3/2)*b*(a + 2*b)**(3/2)/(a + b)**2
  A' = 2*sqrt(a)*b*(2*a + b)/sqrt(a + 2*b)
  A'' = 2*a**(3/2)*b**2*sqrt(a + 2*b)*(2*a + b)**2/(a + b)**4
  A*_O = 2*b*(2*a**3 + 2*a**2*b - 2*a*b**2 + b**3)/(a**(3/2)*sqrt(a + 2*b))
  A'*_O = 2*sqrt(a)*(a + b)**2*(a**2 - 2*a*b + 4*b**2)/(b*(a + 2*b)**(3/2))
  Aped_O = 2*a**(3/2)*b*(a + 2*b)**(3/2)*(a**2 + a*b + b**2)/(a + b)**4
  A'ped_O = sqrt(a)*(a + 2*b)**(3/2)/2

=== k401 at M = O, every ellipse ===
  A'A*_O (I)  = 4*a**2*(a + 2*b)*(a**3 - 2*a**2*b + 2*a*b**2 + 2*b**3)/(b*(2*a + b))    radical-free: True
  A'A*_O (II) = 4*b**2*(2*a + b)*(2*a**3 + 2*a**2*b - 2*a*b**2 + b**3)/(a*(a + 2*b))    radical-free: True
  difference I - II = 4*(a - b)**3*(a + b)**5/(a*b*(a + 2*b)*(2*a + b))
  equals claimed 4(a-b)^3(a+b)^5/(ab(2a+b)(a+2b)): True
  at a=8/5, b=1: 371293/87500  (claim 371293/87500)
  A'A (k106) = 4a^2b^2(2a+b)(a+2b)/(a+b)^2: I ok True, II ok True
  A'/A'' (k113) = (a+b)^4/(ab(2a+b)(a+2b)): I ok True, II ok True
  A/Aped_O (k204 at O) = (a+b)^2/(a^2+ab+b^2): I ok True, II ok True
  A'A'ped_O (k303a at O) = ab(2a+b)(a+2b): I ok True, II ok True

=== foci, a = u^2+1, b = 2u ===
  checks I: {'all_on_ellipse': True, 'closes_after_N': True, 'direction_returns_parallel': True, 'no_earlier_return': True, 'all_sides_tangent': True} 
  checks II: {'all_on_ellipse': True, 'closes_after_N': True, 'direction_returns_parallel': True, 'no_earlier_return': True, 'all_sides_tangent': True}
  orbit I, f1: A'A*_f = -4*(u**2 + 1)**2*(u**2 - 4*u + 1)*(u**2 + 4*u + 1) (rad-free True); equals 4a^2(4b^2-a^2): True; A*_f/A = kappa_6: True
  orbit I, f2: A'A*_f = -4*(u**2 + 1)**2*(u**2 - 4*u + 1)*(u**2 + 4*u + 1) (rad-free True); equals 4a^2(4b^2-a^2): True; A*_f/A = kappa_6: True
  orbit I: k404 A*_f1/Aped_f1 = claimed: True; k204 A/Aped_f1 = claimed: True
  orbit I: A*_f1 == A*_f2: True
  orbit II, f1: A'A*_f = -4*(u**2 + 1)**2*(u**2 - 4*u + 1)*(u**2 + 4*u + 1) (rad-free True); equals 4a^2(4b^2-a^2): True; A*_f/A = kappa_6: True
  orbit II, f2: A'A*_f = -4*(u**2 + 1)**2*(u**2 - 4*u + 1)*(u**2 + 4*u + 1) (rad-free True); equals 4a^2(4b^2-a^2): True; A*_f/A = kappa_6: True
  orbit II: k404 A*_f1/Aped_f1 = claimed: True; k204 A/Aped_f1 = claimed: True
  orbit II: A*_f1 == A*_f2: True
  difference I-II in u: (u - 1)**6*(u + 1)**10/(u*(u**2 + 1)*(u**2 + u + 1)*(u**2 + 4*u + 1)) ; equals (a,b)-formula with a=u^2+1,b=2u: True
