# conicsQuad.py
# Andrew Davison, ad@coe.psu.ac.th, Oct. 2025

"""
This program draws conic sections (ellipse, parabola, or hyperbola) using the general quadratic equation:

A x^2 + Bxy + C y^2 + Dx + Ey + F = 0

The program uses contour plotting to draw smooth curves.
Contour plotting is the standard method for visualizing implicit equations f(x1, x2, x3, ...) = 0

Each conic has both an axis-aligned and a rotated version to illustrate how rotation affects the Bxy term.
"""

import matplotlib.pyplot as plt
import math


def conicEquation(x, y, A, B, C, D, E, F):
  return A * x**2 + B * x * y + C * y**2 + D * x + E * y + F


def classifyConic(A, B, C):
  # Determine conic type from the discriminant
  delta = B**2 - 4 * A * C
  if delta < 0:
    conicType = "Ellipse"
  elif abs(delta) < 1e-10:  # v.close to 0
    conicType = "Parabola"
  else:
    conicType = "Hyperbola"
  print(f"Discriminant delta = {delta:.3f} -> {conicType}")
  return conicType


def plotConic(A, B, C, D, E, F, titleText):
  # Plot a general conic using contour plotting for smooth curves
  conicType = classifyConic(A, B, C)
  xMin, xMax, yMin, yMax = -10, 10, -10, 10
  numPoints = 800  # Higher resolution for smoother curves
  
  # Create grid
  xs = [xMin + i * (xMax - xMin) / (numPoints - 1) for i in range(numPoints)]
  ys = [yMin + j*(yMax-yMin)/(numPoints-1) 
                         for j in range(numPoints)]
  
  # Compute z values for the entire grid
  zs = []
  for y in ys:
    row = []
    for x in xs:
      row.append( conicEquation(x, y, A, B, C, D, E, F))
    zs.append(row)
  
  # Plot the contour where the equation equals zero
  contours = plt.contour(xs, ys, zs, levels=[0], 
                                colors='blue', linewidths=2)
  
  plt.title(titleText)
  plt.xlabel('x')
  plt.ylabel('y')
  plt.axis('equal')
  plt.grid(True, alpha=0.3)
  plt.xlim(xMin, xMax)
  plt.ylim(yMin, yMax)



# --- Conic definitions and plotting ---

print("Select a conic to plot:")
print("1. Ellipse: x^2/a^2 + y^2/b^2 = 1")
print("2. Rotated Ellipse (30 deg)")
print("3. Parabola: y^2 = 4ax")
print("4. Rotated Parabola (45 deg)")
print("5. Hyperbola: x^2/a^2 - y^2/b^2 = 1")
print("6. Rotated Hyperbola (30 deg)")

choice = int(input("Enter your choice (1-6): "))

# Ellipse: x^2/a^2 + y^2/b^2 = 1
if choice == 1:
  a, b = 5, 3
  A = 1 / a**2
  B = 0
  C = 1 / b**2
  D = E = 0
  F = -1
  label = "Ellipse: x^2/a^2 + y^2/b^2 = 1"

# Rotated ellipse
if choice == 2:
  a, b = 5, 3
  angle = math.radians(30)
  A = (math.cos(angle)**2 / a**2) + (math.sin(angle)**2 / b**2)
  B = 2 * math.sin(angle) * math.cos(angle) * (1 / a**2 - 1 / b**2)
  C = (math.sin(angle)**2 / a**2) + (math.cos(angle)**2 / b**2)
  D = E = 0
  F = -1
  label = "Rotated Ellipse (30 deg)"

# Parabola
if choice == 3:
  a = 2
  A = 0
  B = 0
  C = 1
  D = -4 * a
  E = F = 0
  label = "Parabola: y^2 = 4ax"

# Rotated parabola
if choice == 4:
  a = 2
  angle = math.radians(45)
  A = math.sin(angle)**2
  B = -2 * math.sin(angle) * math.cos(angle)
  C = math.cos(angle)**2
  D = -4 * a * (math.cos(angle) - math.sin(angle))
  E = -4 * a * (math.sin(angle) + math.cos(angle))
  F = 0
  label = "Rotated Parabola (45 deg)"

# Hyperbola
if choice == 5:
  a, b = 4, 2
  A = 1 / a**2
  B = 0
  C = -1 / b**2
  D = E = 0
  F = -1
  label = "Hyperbola: x^2/a^2 - y^2/b^2 = 1"

# Rotated hyperbola
if choice == 6:
  a, b = 4, 2
  angle = math.radians(30)
  A = (math.cos(angle)**2 / a**2) - (math.sin(angle)**2 / b**2)
  B = 2 * math.sin(angle) * math.cos(angle) * (1 / a**2 + 1 / b**2)
  C = (math.sin(angle)**2 / a**2) - (math.cos(angle)**2 / b**2)
  D = E = 0
  F = -1
  label = "Rotated Hyperbola (30 deg)"

# Invalid choice
if choice < 1 or choice > 6:
  print("Invalid choice. Please enter a number between 1 and 6.")
  exit()

plt.figure(figsize=(8, 8))
plotConic(A, B, C, D, E, F, label)
plt.show()