# tractrixOrtho.py

'''
Draws circles and finds intersection points between lines and circles.

The tractrix is orthogonal to a set of circles centered on the tractrix's asymptote, all having radius radius R. R is the same as the tractrix's constant length tangent segment.

http://xahlee.info/SpecialPlaneCurves_dir/Tractrix_dir/tractrix.html


the tractrix is a curve perpendicular to a family of equal circles whose centers lie on the X-axis

https://arquimedes.matem.unam.mx/PUEMAC/PUEMAC_2008/rincon/curvas/html/tract.html
'''

import matplotlib.pyplot as plt
from matplotlib.patches import Arc
import math


def lineCircleIntersection(p1, p2, ctr, rad):
  # Direction vector of line
  dx = p2[0] - p1[0]
  dy = p2[1] - p1[1]
  
  # Vector from circle center to first point
  fx = p1[0] - ctr[0]
  fy = p1[1] - ctr[1]
  
  # Quadratic coefficients: a*t^2 + b*t + c = 0
  a = dx * dx + dy * dy
  b = 2 * (fx * dx + fy * dy)
  c = fx * fx + fy * fy - rad * rad
  
  # Calculate discriminant
  disc = b * b - 4 * a * c
  if disc < 0:   # no intersection
    return []
  
  # Tangent (one intersection)
  if disc == 0:
    t = -b / (2 * a)
    x = p1[0] + t * dx
    y = p1[1] + t * dy
    return [(x, y)]
  
  # Two intersections
  sqrtDisc = math.sqrt(disc)
  t1 = (-b - sqrtDisc) / (2 * a)
  t2 = (-b + sqrtDisc) / (2 * a)
  pts = []
  for t in [t1, t2]:
    x = p1[0] + t * dx
    y = p1[1] + t * dy
    pts.append((x, y))
  return pts


def collectPoints(startPt, xs, radius):
  '''
    Start with the point, and define a line using it and the center of the next circle to the right. Calculate the lines intersection with the circle, and store the point. Now repeat using the stored point and the circle to the right of the current one.
  '''
  currPt = startPt
  pts = [currPt]
  for i in range(len(xs) - 1):
    nextCent = (xs[i+1], 0)
    intrs = lineCircleIntersection(currPt, nextCent, nextCent, radius)
    if len(intrs) == 2:
      d1 = (intrs[0][0] - currPt[0])**2 + (intrs[0][1] - currPt[1])**2
      d2 = (intrs[1][0] - currPt[0])**2 + (intrs[1][1] - currPt[1])**2
      currPt = intrs[0] if d1 < d2 else intrs[1]
    elif len(intrs) == 1:
      currPt = intrs[0]
    pts.append(currPt)
  return pts


# --------------------------------------------------------
# centers at y = 0, x from -10 to 10 stepping by 1
xs = list(range(-10, 11, 1))
radius = 5

# Collect points for the top and bottom curves
topPts = collectPoints((xs[0], radius), xs, radius)
bottomPts = collectPoints((xs[0], -radius), xs, radius)

fig, ax = plt.subplots()

# Draw left half-circles
for x in xs:
  halfCircle = Arc((x, 0), 2 * radius, 2 * radius,
                   angle=0, theta1=90, theta2=270)
  ax.add_patch(halfCircle)

# Plot both curves
topXs, topYs = zip(*topPts)
ax.plot(topXs, topYs, 'r-', lw=1)
bottomXs, bottomYs = zip(*bottomPts)
ax.plot(bottomXs, bottomYs, 'b-', lw=1)

ax.set_aspect('equal')
ax.set_xlim(min(xs) - radius - 1, max(xs) + radius + 1)
ax.set_ylim(-radius - 1, radius + 1)
ax.axis('off')

plt.show()