import matplotlib.pyplot as plt
import numpy as np

# 1. 그래프 한글 폰트 및 마이너스 기호 설정
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['axes.unicode_minus'] = False

# ---------------------------------------------------------
# [사용자 입력 부분]
# 직접 입력받고 싶다면 아래 input_data의 주석(#)을 해제하고 사용하세요.
# ---------------------------------------------------------
# print("날짜(t)와 측정값(y)을 공백으로 구분하여 입력하세요 (예: 0, 1, 2, 3...)")
# t_input = list(map(float, input("시간 t (일): ").split()))
# y_input = list(map(float, input("식물 높이 y (cm): ").split()))

# 예시용 샘플 데이터 (사용자가 임의로 수치만 바꾸어 넣으셔도 됩니다)
t_input = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])  # 경과 시간 (일)
y_input = np.array([2.0, 3.5, 7.0, 15.0, 26.0, 37.0, 45.0, 49.0, 50.0])  # 식물 높이 (cm)

t_data = np.array(t_input)
y_data = np.array(y_input)

# 2. 수2 3차 다항함수로 데이터 피팅 (f(t) = a*t^3 + b*t^2 + c*t + d)
# numpy.polyfit을 사용하여 가장 잘 맞는 3차 다항식 계수를 추정합니다.
coefficients = np.polyfit(t_data, y_data, 3)
a, b, c, d = coefficients

# 3. 도함수 f'(t) 계수 구하기 (미분: f'(t) = 3a*t^2 + 2b*t + c)
a_prime = 3 * a
b_prime = 2 * b
c_prime = c

# 4. 연속적인 x축(시간) 구간 생성
t_dense = np.linspace(min(t_data), max(t_data), 500)

# 3차 생장 함수 f(t) 및 2차 도함수 f'(t) 계산
f_t = a * (t_dense**3) + b * (t_dense**2) + c * t_dense + d
f_prime_t = a_prime * (t_dense**2) + b_prime * t_dense + c_prime

# 5. 최대 성장 속도 시점(도함수의 극대점 / 원함수의 변곡점) 수치적 계산
# 2차함수 꼭짓점 축의 방정식: t = -b_prime / (2 * a_prime)
if a_prime != 0:
    t_max_growth = -b_prime / (2 * a_prime)
else:
    t_max_growth = t_dense[np.argmax(f_prime_t)]

# 범위 내에 있는지 확인 후 최대 속도 및 그때의 높이 계산
if min(t_data) <= t_max_growth <= max(t_data):
    max_rate = a_prime * (t_max_growth**2) + b_prime * t_max_growth + c_prime
    height_at_max_rate = (
        a * (t_max_growth**3)
        + b * (t_max_growth**2)
        + c * t_max_growth
        + d
    )
else:
    # 관측 범위 밖일 경우 가장 큰 값 선택
    idx_max = np.argmax(f_prime_t)
    t_max_growth = t_dense[idx_max]
    max_rate = f_prime_t[idx_max]
    height_at_max_rate = f_t[idx_max]

# 결과 출력
print('==================================================')
print(f'추정된 3차 생장 함수: f(t) = {a:.4f}t³ + {b:.4f}t² + {c:.4f}t + {d:.4f}')
print(f'도함수 (성장 속도 함수): f\'(t) = {a_prime:.4f}t² + {b_prime:.4f}t + {c_prime:.4f}')
print('--------------------------------------------------')
print(f'🌱 최대 성장 속도 발생 시점 : t = {t_max_growth:.2f} 일째')
print(f'⚡ 이때의 순간 성장 속도     : {max_rate:.2f} cm/일')
print(f'📏 이때의 식물 높이         : {height_at_max_rate:.2f} cm')
print('==================================================')

# 6. 시각화 (입력 데이터 점 + 모델링 곡선 + 도함수 곡선)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

# [그래프 1] 원함수 f(t) & 실제 측정 데이터
ax1.scatter(
    t_data,
    y_data,
    color='black',
    s=50,
    zorder=5,
    label='Measured Data (입력 데이터)',
)
ax1.plot(
    t_dense,
    f_t,
    color='forestgreen',
    linewidth=2.5,
    label=f'Fitted Model f(t)',
)
ax1.axvline(
    x=t_max_growth,
    color='orange',
    linestyle='--',
    label=f'Max Growth Point (t={t_max_growth:.2f})',
)
ax1.set_ylabel('Plant Height (cm)', fontsize=12)
ax1.set_title(
    'Math II Modeling: User Data Fitting & Growth Rate Analysis', fontsize=14
)
ax1.grid(True, linestyle=':', alpha=0.7)
ax1.legend(loc='upper left')

# [그래프 2] 도함수 f'(t) (순간 성장 속도)
ax2.plot(
    t_dense,
    f_prime_t,
    color='dodgerblue',
    linewidth=2.5,
    label="Growth Rate f'(t)",
)
ax2.scatter(
    [t_max_growth],
    [max_rate],
    color='red',
    s=80,
    zorder=5,
    label=f'Max Rate ({max_rate:.2f} cm/day)',
)
ax2.axvline(x=t_max_growth, color='orange', linestyle='--')
ax2.set_xlabel('Time t (days)', fontsize=12)
ax2.set_ylabel('Growth Rate (cm/day)', fontsize=12)
ax2.grid(True, linestyle=':', alpha=0.7)
ax2.legend(loc='upper right')

plt.tight_layout()
plt.show()