import sympy as sp

# Define symbolic variable
alpha_f = sp.Symbol('alpha_f', real=True)

# Given values from graphs at t = 17.2s
Delta_T_ave = 11.0  # degrees F
d_Delta_T_dt = 1.0  # degrees F/s
lambda_eff = 0.13  # s^-1
dot_lambda_eff = 0.081  # s^-2

# Known parameters
rho_0 = 0.008451  # initial reactivity (845.1 pcm)
alpha_m = -10.6e-5  # in absolute units (-10.6 pcm/F = -10.6e-5 per F)
beta = 0.0065  # delayed neutron fraction

print(f"alpha_m = {alpha_m} per degF = {alpha_m*1e5} pcm/degF")
print()

# Reactivity at t = 17.2s (including BOTH moderator and fuel feedback)
rho = rho_0 + alpha_m * Delta_T_ave + alpha_f * Delta_T_ave

# Rate of change of reactivity
dot_rho = (alpha_m + alpha_f) * d_Delta_T_dt

# Power turning condition (S ≈ 0 at high power):
# dot_rho + (dot_lambda_eff / lambda_eff) * (beta - rho) + lambda_eff * rho = 0

power_turning_eq = dot_rho + (dot_lambda_eff / lambda_eff) * (beta - rho) + lambda_eff * rho

# Solve for alpha_f
solution = sp.solve(power_turning_eq, alpha_f)

print("Power turning equation:")
print(f"dot_rho + (dot_lambda_eff/lambda_eff)*(beta - rho) + lambda_eff*rho = 0")
print()
print("Solving for alpha_f...")
print()

for sol in solution:
    alpha_f_value = float(sol)
    alpha_f_pcm = alpha_f_value * 1e5  # Convert to pcm/degF
    print(f"alpha_f = {alpha_f_value:.6e} per degF")
    print(f"alpha_f = {alpha_f_pcm:.2f} pcm/degF")
    print()

    # Verify the solution
    rho_check = rho_0 + alpha_m * Delta_T_ave + alpha_f_value * Delta_T_ave
    dot_rho_check = (alpha_m + alpha_f_value) * d_Delta_T_dt
    lhs = dot_rho_check + (dot_lambda_eff/lambda_eff)*(beta - rho_check) + lambda_eff*rho_check

    print(f"Verification:")
    print(f"rho at t=17.2s = {rho_check:.6f} ({rho_check*1e5:.1f} pcm)")
    print(f"Power turning equation LHS = {lhs:.6e} (should be ~0)")