
# epicycloid.py
# Andrew Davison, ad@coe.psu.ac.th, Dec. 2025
# https://en.wikipedia.org/wiki/Epicycloid

import math
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

from fractions import Fraction


TIME_STEP = 0.05

def update(frame):
  t = frame * TIME_STEP

  # Center of rolling circle
  centerR = bigR + smallR
  cx = centerR * math.cos(t)
  cy = centerR * math.sin(t)

  # Rolling circle geometry
  rollX, rollY = circlePts(smallR, cx, cy)
  rollingCircleLine.set_data(rollX, rollY)

  # Tracing point
  px, py = epicycloidPt(bigR, smallR, t)
  point.set_data([px], [py])

  pathX.append(px)
  pathY.append(py)
  pathLine.set_data(pathX, pathY)

  return rollingCircleLine, point, pathLine



def epicycloidPt(bigR, smallR, t):
  # Parametric epicycloid equations
  x = (bigR + smallR) * math.cos(t) \
      - smallR * math.cos((bigR + smallR)/smallR * t)
  y = (bigR + smallR) * math.sin(t) \
      - smallR * math.sin((bigR + smallR)/smallR * t)
  return x, y



def circlePts(r, centerX, centerY, steps=200):
  angles = [2 * math.pi * i / steps for i in range(steps + 1)]
  xs = [centerX + r * math.cos(a) for a in angles]
  ys = [centerY + r * math.sin(a) for a in angles]
  return xs, ys


def getRatioInput():
  # Get validated radius ratio from user
  while True:
    try:
      ratio = float(input("Enter R/r radius ratio (>= 1): "))
      if ratio >= 1.0:
        return ratio
      else:
        print("Ratio must be greater than 1.")
    except ValueError:
      print("Please enter a valid number.")


# -----------------------------
bigR = 4.0
ratio = getRatioInput()
smallR = bigR / ratio

ratio_frac = Fraction(ratio).limit_denominator(100)
q = ratio_frac.denominator
# For a hypocycloid with ratio p/q (in lowest terms), the pattern completes after q cycles
frames_needed = int(2 * math.pi * q / TIME_STEP)

fig, ax = plt.subplots()
ax.set_aspect("equal")
ax.set_xlim(-11, 11)
ax.set_ylim(-11, 11)
ax.set_title(f"Epicycloid  R/r={ratio:.2f}")

# Fixed outer circle
bigX, bigY = circlePts(bigR, 0.0, 0.0)
ax.plot(bigX, bigY, color="black", lw=2)

# Rolling circle
rollingCircleLine, = ax.plot([], [], color="blue", lw=2)

# Tracing point
point, = ax.plot([], [], "ro", ms=6)

# Trace path
pathX = []
pathY = []
pathLine, = ax.plot([], [], color="red", lw=1)

anim = FuncAnimation(fig, update, 
           frames=frames_needed, interval=30, blit=True)
plt.show()

