with(plots):
with(plottools):
with(Student[Calculus1]):

# Define the 1-soliton solution
u := (x, t, k, delta) -> (k^2/2)*sech(1/2*(k*x - k^3*t + delta))^2:

# Interactive plot with sliders for t, k, delta
interactive_plot := Explore(
    plot(u(x, t, k, delta), x = -10 .. 10, color = blue, thickness = 2, axes = boxed, labels = ["x", "u(x,t)"]),
    t = -5 .. 5,
    k = 0.5 .. 3,
    delta = -5 .. 5
);

BlaschkeAnalyse := proc(xfun, yfun, tinterval := 0..2*Pi, naam := "Curve", kleur := blue, nulsplits := 3)
local dx, dy, ddx, ddy, kappa, f, i, sol, t1, t2, x1, y1, x2, y2, k1, k2,
      curveplot, tangent_plot, tangents, points, labels, table, solutions, L, norm, vx0, vy0, x0, y0, r, p;

  uses plots, plottools ;
  
  print("�� Analysis of ", naam);

  # Derivatives
  dx := unapply(diff(xfun(t), t), t):
  dy := unapply(diff(yfun(t), t), t):
  ddx := unapply(diff(dx(t), t), t):
  ddy := unapply(diff(dy(t), t), t):

  # Curvature
  kappa := t -> abs(dx(t)*ddy(t) - dy(t)*ddx(t)) / ((dx(t)^2 + dy(t)^2)^(3/2)):

  # Plot the curve
  curveplot := plot([xfun(t), yfun(t), t = tinterval], color = kleur, thickness=2, labels = ["x", "y"]):

  # Difference function f(t) = kappa(t) - kappa(t+Pi)
  f := t -> evalf(kappa(t) - kappa(t + Pi)):

  # Search for roots numerically in split intervals
  solutions := []:
  for i from 0 to nulsplits-1 do
    sol := fsolve(f(t) = 0, t = op(1, tinterval) + i*Pi/nulsplits + 0.01 .. op(1, tinterval) + (i+1)*Pi/nulsplits - 0.01):
    if type(sol, numeric) then
      solutions := [op(solutions), sol]:
    end if;
  end do:

  print("✅ Blaschke pairs found: ", nops(solutions));

  # Visualization elements
  tangents := []:
  points := []:
  labels := []:
  table := []:
  L := 3:

  tangent_plot := proc(t0)
    x0 := evalf(xfun(t0)):
    y0 := evalf(yfun(t0)):
    vx0 := evalf(dx(t0)):
    vy0 := evalf(dy(t0)):
    norm := sqrt(vx0^2 + vy0^2):
    vx0 := vx0 / norm:
    vy0 := vy0 / norm:
    return line([x0 - L*vx0, y0 - L*vy0], [x0 + L*vx0, y0 + L*vy0], color=red, thickness=2):
  end proc:

  for i from 1 to nops(solutions) do
    t1 := evalf(solutions[i]):
    t2 := evalf(t1 + Pi):
    x1 := evalf(xfun(t1)):
    y1 := evalf(yfun(t1)):
    x2 := evalf(xfun(t2)):
    y2 := evalf(yfun(t2)):
    k1 := evalf(kappa(t1)):
    k2 := evalf(kappa(t2)):

    # Tangents
    tangents := [op(tangents), tangent_plot(t1), tangent_plot(t2)]:

    # Points
    points := [op(points),
      pointplot([[x1, y1]], symbol=solidcircle, symbolsize=14, color=red),
      pointplot([[x2, y2]], symbol=solidcircle, symbolsize=14, color=red)
    ]:

    # Labels
    labels := [op(labels),
      textplot([[x1, y1, cat("P", i)]], font=[Times, Bold, 14], align=above, color=red),
      textplot([[x2, y2, cat("P", i, "'")]], font=[Times, Bold, 14], align=above, color=red)
    ]:

    # Table rows
    table := [op(table),
      sprintf("P%d:    t = %.4f   (x,y) = (%.4f, %.4f)   κ = %.4f", i, t1, x1, y1, k1),
      sprintf("P%d':   t = %.4f   (x,y) = (%.4f, %.4f)   κ = %.4f", i, t2, x2, y2, k2),
      "------------------------------------------------------------"
    ]:
  end do:

  # Display table
  print("�� Table of Blaschke pairs:");
  for r in table do
    print(r);
  end do:

  # Show everything
  return display([curveplot, op(tangents), op(points), op(labels)], scaling=constrained,
    title=cat("Blaschke Analysis of ", naam));

end proc:

# Voorbeeldgebruik: superellips
sgn := t -> piecewise(t < 0, -1, t = 0, 0, 1):
a := 3: b := 2: n := 2.5:
x := t -> a * sgn(cos(t)) * abs(cos(t))^n:
y := t -> b * sgn(sin(t)) * abs(sin(t))^n:

BlaschkeAnalyse(x, y, 0..2*Pi, naam="Superellips (n=2.5)", magenta);

BlaschkeAnalyse(
    t -> 5*cos(t),
    t -> (5 + sin(6*t))*sin(t),
    0..2*Pi,
    "Voorbeeldkromme",
    blue,
    3
);