# involuteCycs.py
# Andrew Davison, ad@coe.psu.ac.th, Dec. 2025
'''
Move two points from the crests of adjacent cycloids towards
the cusp located between them. As the points move, draw a tangent
at their current locations which are of the same length as 
the arc-length distance they have covered.

The points at the end of the two tangent segments are collected, and
plotted on the graph, revealing that they form two halves of a
new cycloid.
'''

import math
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


def cycloidPts(ts, a, xOffset):
  xs, ys = zip( *[cycloid(t, a, xOffset) for t in ts])
  return xs, ys


def cycloid(t, a, xOffset=0):
  x = a * (t - math.sin(t)) + xOffset
  y = a * (1 - math.cos(t))
  return (x, y)


def cycloidDf(t, a):
  x = a * (1 - math.cos(t))
  y = a * math.sin(t)
  return (x, y)



def init():
  leftPt.set_data([], [])
  rightPt.set_data([], [])
  leftTangent.set_data([], [])
  rightTangent.set_data([], [])
  leftCurve.set_data([], [])
  rightCurve.set_data([], [])
  return ( leftPt, rightPt,
           leftTangent, rightTangent,
           leftCurve, rightCurve )


def update(frame):
  step = frame * math.pi / steps
  tLeft = math.pi + step
  tRight = math.pi - step

  xLeft, yLeft = cycloid(tLeft, a)
  xRight, yRight = cycloid(tRight, a, xOffset=2.0*math.pi)
  leftPt.set_data([xLeft], [yLeft])
  rightPt.set_data([xRight], [yRight])

  lx, ly, lxEnd, lyEnd = tangentSeg(tLeft, a, 0.0, dir=-1)
  rx, ry, rxEnd, ryEnd = tangentSeg(tRight, a, 2.0*math.pi, dir=1)

  leftTangent.set_data(lx, ly)
  rightTangent.set_data(rx, ry)

  leftTraceX.append(lxEnd); leftTraceY.append(lyEnd)
  rightTraceX.append(rxEnd); rightTraceY.append(ryEnd)

  leftCurve.set_data(leftTraceX, leftTraceY)
  rightCurve.set_data(rightTraceX, rightTraceY)

  return ( leftPt, rightPt,
           leftTangent, rightTangent,
           leftCurve, rightCurve)


def tangentSeg(t, a, xOffset, dir):
  x, y = cycloid(t, a, xOffset)
  dxdt, dydt = cycloidDf(t, a)
  speed = math.hypot(dxdt, dydt)
  if speed == 0.0:
    return [x, x], [y, y], x, y

  ux = dxdt / speed
  uy = dydt / speed
  if dir < 0:
    ux = -ux
    uy = -uy
  arcLength = 4.0 * a * abs(math.cos(t / 2.0))
  xEnd = x + arcLength * ux
  yEnd = y + arcLength * uy
  return [x, xEnd], [y, yEnd], xEnd, yEnd



# -----------------------------------
a = 1.0
tMax = 2.0 * math.pi
steps = 800

ts = [i * tMax / steps for i in range(steps + 1)]
xs1, ys1 = cycloidPts(ts, a, 0.0)
xs2, ys2 = cycloidPts(ts, a, 2.0*math.pi)

fig, ax = plt.subplots()

ax.plot(xs1, ys1)
ax.plot(xs2, ys2)

ax.axhline(0.0)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Cycloid Tangents and Traced Endpoint Curves")
ax.axis("equal")

leftPt, = ax.plot([], [], "ro")
rightPt, = ax.plot([], [], "bo")

leftTangent, = ax.plot([], [], "r")
rightTangent, = ax.plot([], [], "b")

leftCurve, = ax.plot([], [], "r", linewidth=1.5)
rightCurve, = ax.plot([], [], "b", linewidth=1.5)

leftTraceX = []; leftTraceY = []
rightTraceX = []; rightTraceY = []

anim = FuncAnimation(fig, update, frames=steps + 1,
          init_func=init,interval=20, blit=True,repeat=False )

plt.show()