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


import math
import matplotlib.pyplot as plt
from scipy.optimize import brentq


def catenoidArea(R, H):
  # Solve for catenoid parameter a using brentq()
  z = H/2
  aMin = R * math.exp(-z/R)
  a = brentq(catenoidEq, aMin, R, args=(R, z))

  return math.pi * a * (H + a * math.sinh(H / a))


def catenoidEq(a, R, z):
  return a * math.cosh(z/a) - R


def cylinderArea(R, H):
  return 2*math.pi * R * H


# ---------------------------------------------------------
R = 1.0      # fixed radius
CRITICAL_HR = 1.32548 
# Catenoid exists only if H/R < CRITICAL_HR

Hmin = 0.1   # avoid H = 0
Hmax = CRITICAL_HR*R
numSteps = 80

hs = [(Hmin + (Hmax-Hmin)*(i/numSteps)) for i in range(numSteps)]
catAreas = [ catenoidArea(R, h) for h in hs]
cylAreas = [ cylinderArea(R, h) for h in hs]

plt.figure()

plt.plot(hs, catAreas, label="Catenoid Area")
plt.plot(hs, cylAreas, label="Cylinder Area")

plt.xlabel("H (height)")
plt.ylabel("Surface Area")
plt.title("Surface Area vs H for Catenoid and Cylinder")
plt.legend()
plt.grid(True)

plt.show()