import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# ----------------------------
# Physical constants
# ----------------------------
g = 9.81
a_out = 0.01
H0 = 1.0
V_target = 1.0
# ----------------------------
def empty_time(A_func):
    integrand = lambda h: A_func(h) / (a_out * np.sqrt(2*g*h))
    return quad(integrand, 1e-8, H0)[0]

# ----------------------------
# Sweep r values
# ----------------------------
r_values = np.linspace(1, 10, 200)
ratios = []

for r in r_values:

    # Area scaling to keep total volume constant
    A_small = 3*V_target / (H0 * (1 + r + np.sqrt(r)))
    A_large = r * A_small

    # Tank A: small bottom → large top
    def A_A(h):
        return A_small + (A_large - A_small) * h / H0

    # Tank B: large bottom → small top
    def A_B(h):
        return A_large - (A_large - A_small) * h / H0

    T_A = empty_time(A_A)
    T_B = empty_time(A_B)

    ratios.append(T_B / T_A)

# ----------------------------
# Plot
# ----------------------------
plt.figure(figsize=(8,5))
plt.plot(r_values, ratios, color='blue')
plt.axhline(1, linestyle='--', color='black')
plt.xlabel("$r={(A_{top}/ A_{bottom})}_{tankA}$")
plt.ylabel("$T_B / T_A$")
plt.title("Emptying Time Ratio vs Cross-section Ratio")
plt.grid(True)
plt.show()