import numpy as np

# Parámetros del problema
G_tot = 1000  # Irradiación solar total (W/m^2)
A_conc = 10  # Área del concentrador (m^2)
alpha_conc = 0.8  # Absorbancia del concentrador
sigma = 5.67e-8  # Constante de Stefan-Boltzmann (W/m^2 K^4)
epsilon_pv = 0.9  # Emisividad de las células fotovoltaicas
A_pv = 5  # Área de las células fotovoltaicas (m^2)
h_conc = 10  # Coeficiente de transferencia de calor por convección (W/m^2 K)
epsilon_conc = 0.85  # Emisividad del concentrador
Tem_sky = 300  # Temperatura del cielo (K)
Tem_a = 298  # Temperatura ambiente (K)

# Funciones para calcular las tasas de transferencia de calor
def q_solar_conc(Tem_conc):
    return G_tot * A_conc * alpha_conc

def q_rad_pv(Tem_conc):
    return sigma * epsilon_pv * A_pv * (Tem_conc**4 - Tem_sky**4)

def q_conv_conc(Tem_conc):
    return 2 * A_conc * h_conc * (Tem_conc - Tem_a)

def q_rad_conc(Tem_conc):
    return sigma * epsilon_conc * A_conc * (Tem_conc**4 - Tem_sky**4)

# Método de Runge-Kutta de cuarto orden
def runge_kutta_4(Tem_conc, dt):
    k1 = (q_solar_conc(Tem_conc) + q_rad_pv(Tem_conc)) - (q_conv_conc(Tem_conc) + q_rad_conc(Tem_conc))
    k2 = (q_solar_conc(Tem_conc + k1*dt/2) + q_rad_pv(Tem_conc + k1*dt/2)) - (q_conv_conc(Tem_conc + k1*dt/2) + q_rad_conc(Tem_conc + k1*dt/2))
    k3 = (q_solar_conc(Tem_conc + k2*dt/2) + q_rad_pv(Tem_conc + k2*dt/2)) - (q_conv_conc(Tem_conc + k2*dt/2) + q_rad_conc(Tem_conc + k2*dt/2))
    k4 = (q_solar_conc(Tem_conc + k3*dt) + q_rad_pv(Tem_conc + k3*dt)) - (q_conv_conc(Tem_conc + k3*dt) + q_rad_conc(Tem_conc + k3*dt))
    
    Tem_conc_next = Tem_conc + (k1 + 2*k2 + 2*k3 + k4) * dt / 6
    
    return Tem_conc_next

# Condiciones iniciales y parámetros del problema
Tem_conc_inicial = 350  # Temperatura inicial del concentrador (K)
dt = 0.01  # Tamaño del paso en el tiempo
tiempo_final = 10.0  # Tiempo final de simulación

# Iterar en el tiempo para resolver la ecuación para Tem_conc
Tem_conc = Tem_conc_inicial
tiempo_actual = 0
while tiempo_actual < tiempo_final:
    Tem_conc = runge_kutta_4(Tem_conc, dt)
    tiempo_actual += dt

# En este punto, Tem_conc contiene la temperatura estimada del concentrador solar en el tiempo final deseado
print("Temperatura final del concentrador solar:", Tem_conc, "K")