# rollParabola.py
# Andrew Davison, Nov 2025, ad@coe.psu.ac.th

'''
Animates a parabola rolling along the x-axis. The parabola rotates clockwise
when moving right and counterclockwise when moving left. It reverses direction
when the tip of either arm touches the x-axis. The focus point is tracked and
its path is displayed. Uses pre-calculated arc lengths for proper rolling motion.

In addition, the curve y = 1/(4*a) * cosh(4*a*x) is drawn, and visually we can see that the focus path follows this curve exactly.

The parabola's shape is controlled by the equation y = ax^2, where the parameter 'a' determines its characteristics:

The coefficient 'a':
- Larger values (like a = 1.0): Make the parabola narrower and steeper. The sides rise more quickly, creating a tall, thin U-shape
- Smaller values (like a = 0.2): Make the parabola wider and flatter. The sides rise more gradually, creating a shallow, broad U-shape
- Very small values (like a = 0.05): Create an extremely wide parabola that's almost flat near the vertex

The range parameter 't' (in the code: tMin to tMax):
- Controls how far along the parabola's arms we draw
- Larger range like [-6, 6] extends the arms further out horizontally
- Smaller range like [-2, 2] only shows the parabola near its vertex
- Since y = ax^2, the actual height at the endpoints is a × t^2

The focal length (1/(4a)):
- Gets smaller as 'a' increases
- For a = 1.0, the focus is only 0.25 units above the vertex
- For a = 0.2, the focus is 1.25 units above the vertex
- This is why narrower parabolas have their focus closer to the vertex

The focal length formula 1/(4a) comes from the mathematical definition of a parabola as the set of all points equidistant from a fixed point (the focus) and a fixed line (the directrix).

Derivation:

For a parabola with equation y = ax^2 and vertex at the origin:

1. The focus is located at point (0, f) where f is the focal length we want to find
2. The directrix is a horizontal line at y = -f (same distance below the vertex as the focus is above)
3. By definition, any point (x, y) on the parabola must be equidistant from the focus and the directrix

Setting up the equation:
- Distance from point (x, y) to focus (0, f): √(x^2 + (y - f)^2)
- Distance from point (x, y) to directrix y = -f: |y - (-f)| = y + f

Setting them equal:
√(x^2 + (y - f)^2) = y + f

Squaring both sides:
x^2 + (y - f)^2 = (y + f)^2

Expanding:
x^2 + y^2 - 2yf + f^2 = y^2 + 2yf + f^2

Simplifying:
x^2 = 4yf

Solving for y:
y = x^2/(4f)

Comparing with y = ax^2:
Since our parabola is y = ax^2, we can see that:
- a = 1/(4f)
- Therefore: f = 1/(4a)

This is why the focal length is always 1/(4a) for a parabola of the form y = ax^2.

'''

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


TOL = 0.2

def catenary(x, a):
  return 1/(4*a) * math.cosh(4*a*x)


def calcLengths(xs, ys):
  # Calculate arc lengths from vertex to each point
  # Vertex is at center of symmetric range
  nPts = len(xs)
  vertIdx = nPts//2
  arcLens = [0.0] * nPts
  
  # Calculate arc lengths to the right of vertex (positive)
  for i in range(vertIdx + 1, nPts):
    dx = xs[i] - xs[i-1]
    dy = ys[i] - ys[i-1]
    segLen = math.sqrt(dx*dx + dy*dy)
    arcLens[i] = arcLens[i-1] + segLen
  
  # Calculate arc lengths to the left of vertex (negative)
  for i in range(vertIdx - 1, -1, -1):
    dx = xs[i] - xs[i+1]
    dy = ys[i] - ys[i+1]
    segLen = math.sqrt(dx*dx + dy*dy)
    arcLens[i] = arcLens[i+1] - segLen
  return arcLens



def rotatePt(x, y, angle):
  cosA = math.cos(angle)
  sinA = math.sin(angle)
  xRot = x*cosA - y*sinA
  yRot = x*sinA + y*cosA
  return xRot, yRot


def init():
  parLine.set_data([], [])
  focusPt.set_data([], [])
  focusPath.set_data([], [])
  contactPt.set_data([], [])
  return parLine, focusPt, focusPath, contactPt


def animate(frame):
  global rotAngle, direction
  
  rotAngle += direction * angleStep
  
  # Rotate each point
  parXRot = []
  parYRot = []
  for i in range(len(parLocalX)):
    xr, yr = rotatePt(parLocalX[i], parLocalY[i], rotAngle)
    parXRot.append(xr)
    parYRot.append(yr)
  
  # Find the lowest point after rotation (contact point)
  minYIdx = 0
  minYVal = parYRot[0]
  for i in range(1, len(parYRot)):
    if parYRot[i] < minYVal:
      minYVal = parYRot[i]
      minYIdx = i
  
  '''
   Position contact point at its arc length distance from origin
   Arc length is already signed: negative for left branch, 
   positive for right branch
  '''
  arcLenAtContact = arcLensVert[minYIdx]
  currX = arcLenAtContact
  
  contactXLocal = parXRot[minYIdx]
  contactYLocal = parYRot[minYIdx]
  
  # Calculate offsets
  xOffset = currX - contactXLocal
  yOffset = -contactYLocal
  
  # translate
  parX = [x + xOffset for x in parXRot]
  parY = [y + yOffset for y in parYRot]
  
  # Check if the tips touch the x-axis
  # in order to change direction
  leftEndY = parY[0]
  rightEndY = parY[-1]
  if leftEndY < TOL or rightEndY < TOL:
    direction *= -1
  
  # move focus position
  focusXRot, focusYRot = rotatePt(0, focalHeight, rotAngle)
  focusX = focusXRot + xOffset
  focusY = focusYRot + yOffset
  focusXs.append(focusX)
  focusYs.append(focusY)
  
  # store contact point
  contactX = currX
  contactY = 0
  
  # Update plot elements
  parLine.set_data(parX, parY)
  focusPt.set_data([focusX], [focusY])
  focusPath.set_data(focusXs, focusYs)
  contactPt.set_data([contactX], [contactY])
  
  return parLine, focusPt, focusPath, contactPt


# ----------------------------------------
# Parabola parameters
a = 0.4
focalHeight = 1/(4*a)

# Create parabola curve in local coords
numPts = 200
tMin, tMax = -14, 14
ts = [tMin + (tMax-tMin)*i /(numPts-1) for i in range(numPts)]
parLocalX = ts
parLocalY = [a * t**2 for t in ts]

cats = [catenary(t,a) for t in ts]

arcLensVert = calcLengths(parLocalX, parLocalY)

# Setup figure and axis
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_xlim(-15, 15)
ax.set_ylim(-1, 20)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='k', lw=0.5)
ax.axvline(x=0, color='k', lw=0.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Rolling Parabola with Focus Path')

ax.plot(ts, cats, 'k--', lw=1, alpha=0.6, label='Catenary')

# Initialize plot elements
parLine, = ax.plot([], [], 'b-', lw=2, label='Parabola')
focusPt, = ax.plot([], [], 'ro', markersize=8, label='Focus')
focusPath, = ax.plot([], [], 'r', lw=1, alpha=0.6, label='Focus Path')
contactPt, = ax.plot([], [], 'go', markersize=6, label='Contact Pt')

# Storage for focus path
focusXs = []
focusYs = []

# Animation state
rotAngle = 0.0
direction = 1   # turn left initially
angleStep = 0.02

anim = FuncAnimation(fig, animate, init_func=init, frames=2000, 
                     interval=30, blit=True, repeat=True)

ax.legend(loc='upper right')
plt.tight_layout()
plt.show()