import numpy as np
from scipy.optimize import fsolve

# Parámetros
G_tot = 1000  # Irradiancia total en W/m²
G_b = 800  # Irradiancia directa en W/m²
A_conc = 13.77  # Área del concentrador en m²
A_pv = 1.2  # Área del PV en m²
A_abs = 0.6  # Área del absorbedor en m²
alpha_conc = 0.03  # Absorptividad del concentrador
alpha_pv = 0.97  # Absorptividad del PV
alpha_abs = 0.9  # Absorptividad del absorbedor
sigma = 5.67e-8  # Constante de Stefan-Boltzmann en W/(m²·K⁴)
epsilon_pv = 0.2  # Emisividad del PV
epsilon_abs = 0.2  # Emisividad del absorbedor
epsilon_conc = 0.3  # Emisividad del concentrador
eta_opt = 0.9  # Eficiencia óptica
IAM_th = 1.0  # Modificador angular térmico
IAM_elec = 0.95  # Modificador angular eléctrico
eta_pv = 0.15  # Eficiencia del PV
h_conc = 10  # Coeficiente de convección en W/(m²·K)
R_conv_pv = 0.1  # Resistencia convectiva del PV
R_conv_abs = 0.1  # Resistencia convectiva del absorbedor
R_cond_pv = 0.05  # Resistencia conductiva del PV
R_cond_abs = 0.05  # Resistencia conductiva del absorbedor
R_cond_sub = 0.02  # Resistencia conductiva del substrato
C_p_f = 4186  # Capacidad calorífica del fluido en J/(kg·K)
m_dot_f = 0.15  # Flujo másico del fluido en kg/s
Tem_a = 298  # Temperatura ambiente en K
Tem_sky = 298  # Temperatura del cielo en K
Tem_f_in = 343  # Temperatura de entrada del fluido en K

# Definir las ecuaciones internas
def equations(vars):
    Tem_conc, Tem_pv, Tem_abs, Tem_f_out, Tem_sub = vars
    
    # Ecuaciones de balance de energía
    q_solar_conc = G_tot * A_conc * alpha_conc
    q_rad_pv = sigma * epsilon_pv * A_pv * (Tem_pv**4 - Tem_conc**4)
    q_conv_conc_front = A_conc * h_conc * (Tem_conc - Tem_a)
    q_conv_conc_back = A_conc * h_conc * (Tem_conc - Tem_a)
    q_rad_conc = sigma * epsilon_conc * A_conc * (Tem_conc**4 - Tem_sky**4)
    
    # Ecuaciones de componentes individuales
    eq1 = q_solar_conc + q_rad_pv - (q_conv_conc_front + q_conv_conc_back + q_rad_conc)
    
    q_solar_pv = G_b * A_pv * eta_opt * IAM_th * alpha_pv
    q_solar_abs = G_tot * A_abs * alpha_abs
    P_ele_pv = G_b * A_pv * eta_opt * IAM_elec * eta_pv
    q_rad_abs = sigma * epsilon_abs * A_abs * (Tem_abs**4 - Tem_sky**4)
    q_conv_pv = (Tem_pv - Tem_a) * R_conv_pv
    q_conv_abs = (Tem_abs - Tem_a) * R_conv_abs
    q_cond_pv_x_sub = (Tem_pv - Tem_sub) / (R_cond_pv + R_cond_sub)
    q_cond_abs_x_sub = (Tem_abs - Tem_sub) / (R_cond_abs + R_cond_sub)
    q_f = m_dot_f * C_p_f * (Tem_f_out - Tem_f_in)
    
    # Ecuaciones de balance de energía para el PV y el absorbedor
    eq2 = q_solar_pv - P_ele_pv - q_conv_pv - q_rad_pv + q_solar_abs - q_conv_abs - q_rad_abs - q_f
    eq3 = q_solar_abs - q_conv_abs - q_rad_abs - q_cond_abs_x_sub
    eq4 = m_dot_f * C_p_f * (Tem_f_out - Tem_f_in) - q_cond_pv_x_sub - q_cond_abs_x_sub
    
    # Condición de equilibrio térmico para el fluido
    eq5 = q_cond_pv_x_sub + q_cond_abs_x_sub - q_f
    
    return [eq1, eq2, eq3, eq4, eq5]

# Estimación inicial para las variables
initial_guess = [350, 350, 350, 310, 300]  # Estimaciones iniciales para [Tem_conc, Tem_pv, Tem_abs, Tem_f_out, Tem_sub]

# Resolver el sistema de ecuaciones no lineales
solution = fsolve(equations, initial_guess)

Tem_conc_solution, Tem_pv_solution, Tem_abs_solution, Tem_f_out_solution, Tem_sub_solution = solution

# Convertir las soluciones a grados Celsius
Tem_conc_solution_C = Tem_conc_solution - 273.15
Tem_pv_solution_C = Tem_pv_solution - 273.15
Tem_abs_solution_C = Tem_abs_solution - 273.15
Tem_f_out_solution_C = Tem_f_out_solution - 273.15
Tem_sub_solution_C = Tem_sub_solution - 273.15

print(f"La temperatura de equilibrio del concentrador es: {Tem_conc_solution:.2f} K ({Tem_conc_solution_C:.2f} °C)")
print(f"La temperatura de equilibrio del PV es: {Tem_pv_solution:.2f} K ({Tem_pv_solution_C:.2f} °C)")
print(f"La temperatura de equilibrio del absorbedor es: {Tem_abs_solution:.2f} K ({Tem_abs_solution_C:.2f} °C)")
print(f"La temperatura de salida del fluido es: {Tem_f_out_solution:.2f} K ({Tem_f_out_solution_C:.2f} °C)")
print(f"La temperatura del substrato es: {Tem_sub_solution:.2f} K ({Tem_sub_solution_C:.2f} °C)")