# confocalPara.py
# Andrew Davison, ad@coe.psu.ac.th, Oct. 2025
'''
Generate confocal parabolas with the same focus 
(f,0) but different directrices:
* some open right (d < f)
* some open left (d > f)
'''

import matplotlib.pyplot as plt
import bisect  # for sorted list insertion


def parabolaPlot(xMin, xMax, f, d):
  '''
  Returns lists of x and y values for a parabola with focus (f,0) and vertical directrix x = d, between xMin and xMax. The parabola equation is y^2 = 4p(x - (f - p)), 
  where p = (f - d)/2.
  '''
  p = (f - d) / 2
  vX = f - p   # vertex of parabola (where crosses y-axis)

  if xMax == xMin:
    return [], []

  stepSize = (xMax - xMin) / steps
  xsRange = [xMin + i * stepSize for i in range(steps+1)]
  # Insert vX in sorted order if within range;
  # avoids gaps in the plotted points
  if xMin <= vX <= xMax:
    bisect.insort(xsRange, vX)

  xs =  []; ys = []
  for x in xsRange:
    ySqu = 4 * p * (x - vX)
    if ySqu >= 0:  # only store positive y's
      ys.append(ySqu ** 0.5)
      xs.append(x)
  return xs, ys


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

f = 3       # focus at (3,0)
ds = [-1, 1, 2.5, 4, 6, 8] # directrix positions
# d < f --> opens right, d > f --> opens left
steps = 100

plt.figure(figsize=(5, 8))

for d in ds:
  xs, ys = parabolaPlot(-5, 10, f, d)
  plt.plot(xs, ys, 'b')
  plt.plot(xs, [-y for y in ys], 'b')
  if d < f:  # opens right
    plt.text(xs[-1]-1.5, ys[-1]+0.5, 'd='+str(d)) 
  else:   # opens left
    plt.text(xs[0], ys[0]+0.5, 'd='+str(d)) 

plt.plot(f, 0, 'ro')  # Focal point

plt.axhline(0, color='k', ls='-.', lw=1)
plt.axvline(0, color='k', ls='-.', lw=1)
plt.axis('equal')
# plt.axis('off')
plt.title("Confocal Parabolas Facing Left and Right")
plt.show()