program Newton2Vars_HTF
    implicit none
    real(8) :: T_abs, T_sub
    real(8) :: T_a, T_sky, T_PV, T_in_HTF, T_PTC
    real(8) :: G_GHI, A_abs, alpha_abs, epsilon_abs
    real(8) :: sigma, R_conv_abs, R_cond_abs, R_cond_sub, R_cond_PV
    real(8) :: m_dot_HTF, C_p_HTF, h_HTF
    real(8) :: f1, f2, J11, J12, J21, J22
    real(8) :: dT_abs, dT_sub, detJ
    real(8) :: tol, UA_HTF, NTU, eff_HTF
    integer :: iter, max_iter

    sigma = 5.67e-8

    ! ==== Datos conocidos ====
    G_GHI      = 1000.0d0
    A_abs      = 0.6d0
    alpha_abs  = 0.9d0
    epsilon_abs= 0.2d0

    T_a        = 298.2d0
    T_sky      = 298.2d0
    T_PV       = 355.5d0
    T_in_HTF   = 343.2d0
    T_PTC      = 301.5d0

    R_conv_abs = 0.0464d0
    R_cond_abs = 0.00002439d0
    R_cond_sub = 0.00004623d0
    R_cond_PV  = 0.00005d0  ! <<--- asignar valor apropiado

    m_dot_HTF  = 0.15d0
    C_p_HTF    = 4187.0d0
    h_HTF      = 1665.0d0

    ! ==== Configuración Newton–Raphson ====
    T_abs = 330.0d0
    T_sub = 320.0d0
    tol = 1.0d-6
    max_iter = 50

    do iter = 1, max_iter

        ! ---- Calcular UA y efectividad del HTF ----
        UA_HTF = 1.0d0 / (1.0d0/h_HTF + R_cond_sub)
        NTU = (UA_HTF * (3.1416d0 * 0.03d0 * 10.0d0)) / (m_dot_HTF * C_p_HTF)
        eff_HTF = 1.0d0 - exp(-NTU)

        ! ---- Funciones ----
        f1 = G_GHI*A_abs*alpha_abs &
             - (T_abs - T_a)/R_conv_abs &
             - epsilon_abs*sigma*A_abs*(T_abs**4 - T_sky**4) &
             - (T_abs - T_sub)/(R_cond_abs + R_cond_sub)

        f2 = (T_PV - T_sub)/(R_cond_PV + R_cond_sub) &
             + (T_abs - T_sub)/(R_cond_abs + R_cond_sub) &
             - eff_HTF * m_dot_HTF * C_p_HTF * (T_sub - T_in_HTF)

        ! ---- Jacobiano ----
        ! Parciales de f1
        J11 = -1.0d0/R_conv_abs - 4.0d0*epsilon_abs*sigma*A_abs*T_abs**3 - 1.0d0/(R_cond_abs + R_cond_sub)
        J12 =  1.0d0/(R_cond_abs + R_cond_sub)

        ! Parciales de f2
        J21 =  1.0d0/(R_cond_abs + R_cond_sub)
        J22 = -1.0d0/(R_cond_PV + R_cond_sub) - 1.0d0/(R_cond_abs + R_cond_sub) - eff_HTF*m_dot_HTF*C_p_HTF

        ! ---- Determinante ----
        detJ = J11*J22 - J12*J21

        ! ---- Correcciones ----
        dT_abs = (-f1*J22 + f2*J12) / detJ
        dT_sub = (-J11*f2 + J21*f1) / detJ

        ! ---- Actualizar ----
        T_abs = T_abs + dT_abs
        T_sub = T_sub + dT_sub

        ! ---- Convergencia ----
        if (max(abs(dT_abs), abs(dT_sub)) < tol) exit

    end do

    ! ==== Resultados ====
    print *, "Iteraciones: ", iter
    print *, "T_abs (K): ", T_abs, "  T_abs (°C): ", T_abs - 273.15d0
    print *, "T_sub (K): ", T_sub, "  T_sub (°C): ", T_sub - 273.15d0

end program Newton2Vars_HTF