# cycloidA.py
# Andrew Davison, ad@coe.psu.ac.th, Dec. 2025
"""
Cycloid Visualization with slider to adjust the 
rolling circle radius 'a'.

Shows:
  * Area under one arch (between two consecutive cusps):
    - Analytical formula: 3*pi*a^2
    - Numerical verification using Simpson's rule integration
  
  * Arc length of one arch:
    - Analytical formula: 8a
    - Numerical verification using Simpson's rule integration
    - Uses the arc length formula: L = integral sqrt((dx/dt)^2 + (dy/dt)^2) dt
"""

import math
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from scipy.integrate import simpson


MAX_X = 27
MAX_Y = 5
nPts = 2000


def cycloid(ts, a=1):
  xs = [a * (t - math.sin(t)) for t in ts]
  ys = [a * (1 - math.cos(t)) for t in ts]
  return xs, ys


def cycloidDf(t, a=1):
  # Derivative of the cycloid
  dx_dt = a * (1 - math.cos(t))
  dy_dt = a * math.sin(t)
  return (dx_dt, dy_dt)


def calcArea(a):
  # Calculate area under one arch (between first two cusps)
  analyticalArea = 3 * math.pi * a * a
  
  # Numerical integration using Simpson's rule
  # For one arch, t goes from 0 to 2 pi
  ts = [i * (2*math.pi)/(nPts-1) for i in range(nPts)]
  
  # Get x and y values from the cycloid function
  xs, ys = cycloid(ts, a)
  
  # Use Simpson's rule to integrate y with respect to x
  # We need dx values for integration
  dxVals = [xs[i+1] - xs[i] for i in range(len(xs)-1)]
  
  # For Simpson's rule with parametric curves, integrate y * (dx/dt) * dt
  # dx/dt = a(1 - cos(t))
  dxdt = [a * (1 - math.cos(t)) for t in ts]
  integrand = [ys[i] * dxdt[i] for i in range(len(ys))]
  dt = (2 * math.pi) / (nPts - 1)
  numericalArea = simpson(integrand, dx=dt)
  
  return analyticalArea, numericalArea


def calcArcLen(a):
  # Calculate arc length of one arch (between first two cusps)
  analyticalLen = 8 * a
  
  # Numerical integration using Simpson's rule
  # For one arch, t goes from 0 to 2 pi
  ts = [i * (2*math.pi)/(nPts-1) for i in range(nPts)]
  
  # Arc length formula: L = integral of sqrt((dx/dt)^2 + (dy/dt)^2) dt
  # Use cycloidDf to get derivatives
  integrand = []
  for t in ts:
    dx_dt, dy_dt = cycloidDf(t, a)
    integrand.append(math.sqrt(dx_dt**2 + dy_dt**2))
  
  dt = (2 * math.pi) / (nPts - 1)
  numericalLen = simpson(integrand, dx=dt)
  
  return analyticalLen, numericalLen



def update(val):
  drawCanvas(a_slider.val)
  fig.canvas.draw_idle()


def drawCanvas(a):
  ax.clear()

  xs, ys = cycloid(ts, a)
  ax.plot(xs, ys, lw=2)

  ax.set_title("Cycloid with Adjustable Parameter a")
  ax.set_xlabel("x")
  ax.set_ylabel("y")
  ax.set_xlim(0, MAX_X)
  ax.set_ylim(0, MAX_Y)

  # X-axis integer ticks
  xticks = list(range(MAX_X + 1))
  ax.set_xticks(xticks)
  ax.set_xticklabels([str(i) for i in xticks])

  # X-axis pi ticks (minor)
  pi_xticks = [k * math.pi for k in range(9)]
  ax.set_xticks(pi_xticks, minor=True)

  # Make pi ticks longer and thicker
  ax.tick_params( axis="x", which="minor", length=12, width=2)
  for k, x in enumerate(pi_xticks):
    label = "0" if k == 0 else f"{k}π"
    ax.text(x, -0.15, label,
      transform=ax.get_xaxis_transform(),
      ha="center", va="top")
    ax.axvline(x=x, color="red",linestyle="--", 
               lw=0.8, alpha=0.6, zorder=0)

  # Y-axis integer ticks
  yticks = list(range(MAX_Y + 1))
  ax.set_yticks(yticks)
  ax.set_yticklabels([str(ytick) for ytick in yticks])
  for ytick in yticks:      # Horizontal lines
    ax.axhline(y=ytick, color="red",linestyle="--", 
               lw=0.8, alpha=0.6, zorder=0)

  # Calculate and display the area and arc len
  analyticalArea, numericalArea = calcArea(a)
  second_cusp_x = 2 * math.pi * a
  info = f"Area (0 to {second_cusp_x:.3f}):\n"
  info += f"Analytical: {analyticalArea:.3f} (3πa²)\n"
  info += f"Numerical: {numericalArea:.3f} (Simpson)\n\n"
  
  analyticalLen, numericalLen = calcArcLen(a)
  info += f"Arc Len (one arch):\n"
  info += f"Analytical: {analyticalLen:.3f} (8a)\n"
  info += f"Numerical: {numericalLen:.3f} (Simpson)"
  
  ax.text(0.98, 0.97, info,
          transform=ax.transAxes,
          verticalalignment='top',
          horizontalalignment='right')



# ----------------------------------

step = (8 * math.pi) / 1999
ts = [i * step for i in range(2000)]

a0 = 1.0

fig, ax = plt.subplots(figsize=(11, 5))
plt.subplots_adjust(bottom=0.35)

drawCanvas(a0)

slider_ax = fig.add_axes([0.15, 0.1, 0.7, 0.05])
a_slider = Slider(ax=slider_ax, label="a", 
                  valmin=0.2, valmax=2.0,valinit=a0, valstep=0.05)
a_slider.on_changed(update)

plt.show()