program Newton3Vars_HTF
    implicit none
    integer, parameter :: dp = selected_real_kind(15,300)
    ! variables incógnitas
    real(dp) :: T_abs, T_sub, T_out
    ! datos y constantes
    real(dp) :: T_a, T_sky, T_PV, T_in_HTF, T_PTC
    real(dp) :: q_dot_sol_PV, P_PV
    real(dp) :: G_GHI, A_abs, alpha_abs, epsilon_abs, A_PV, epsilon_PV
    real(dp) :: sigma, R_conv_abs, R_cond_abs, R_cond_sub, R_cond_PV, R_conv_PV
    real(dp) :: m_dot_HTF, C_p_HTF, h_HTF
    real(dp) :: UA, A_pipe, NTU, eff_HTF
    ! funciones y jacobiano
    real(dp) :: f1, f2, f3
    real(dp) :: J(3,3)
    real(dp) :: rhs(3), dx(3)
    real(dp) :: tol
    integer :: iter, max_iter, i

    ! Inicializar datos (tus valores)
    sigma = 5.67d-8
    G_GHI      = 1000.0d0
    q_dot_sol_PV = 540.0d0
    P_PV = 1722.0d0

    A_abs      = 0.6d0
    alpha_abs  = 0.9d0
    epsilon_abs= 0.2d0
    A_PV       = 1.2d0
    epsilon_PV = 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
    R_conv_PV  = 0.03281d0

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

    ! geometría tubería
    A_pipe = 3.1416d0 * 0.03d0 * 10.0d0

    ! calcular UA, NTU, eff (constantes)
    UA = 1.0d0 / (1.0d0/h_HTF + R_cond_sub)
    NTU = (UA * A_pipe) / (m_dot_HTF * C_p_HTF)
    eff_HTF = 1.0d0 - exp(-NTU)

    ! iniciales (elige valores razonables)
    T_abs = 330.0d0
    T_sub = 320.0d0
    T_out = 345.0d0

    tol = 1.0d-8
    max_iter = 100

    do iter = 1, max_iter
        ! --- calcular f ---
        f1 = q_dot_sol_PV - P_PV - (T_PV - T_a)/R_conv_PV &
             - epsilon_PV * sigma * A_PV * (T_PV**4 - T_PTC**4) &
             + 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) &
             - eff_HTF * m_dot_HTF * C_p_HTF * (T_sub - T_in_HTF)

        f2 = 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)

        f3 = m_dot_HTF * C_p_HTF * (T_out - T_in_HTF) &
             - eff_HTF * m_dot_HTF * C_p_HTF * (T_sub - T_in_HTF)

        ! --- montar J (parciales) ---
        J(1,1) = -1.0d0/R_conv_abs - 4.0d0*epsilon_abs*sigma*A_abs*T_abs**3 - 1.0d0/(R_cond_abs + R_cond_sub)
        J(1,2) =  1.0d0/(R_cond_abs + R_cond_sub) - eff_HTF * m_dot_HTF * C_p_HTF
        J(1,3) =  0.0d0

        J(2,1) = -1.0d0/R_conv_abs - 4.0d0*epsilon_abs*sigma*A_abs*T_abs**3 - 1.0d0/(R_cond_abs + R_cond_sub)
        J(2,2) =  1.0d0/(R_cond_abs + R_cond_sub)
        J(2,3) =  0.0d0

        J(3,1) =  0.0d0
        J(3,2) = - eff_HTF * m_dot_HTF * C_p_HTF
        J(3,3) =   m_dot_HTF * C_p_HTF

        ! --- rhs = -f ---
        rhs(1) = -f1
        rhs(2) = -f2
        rhs(3) = -f3

        ! --- resolver J * dx = rhs por eliminación gaussiana (3x3) ---
        call solve3x3(J, rhs, dx)

        ! actualizar
        T_abs = T_abs + dx(1)
        T_sub = T_sub + dx(2)
        T_out = T_out + dx(3)

        ! condición de paro
        if (max(abs(dx(1)), max(abs(dx(2)), abs(dx(3)))) < tol) exit

        ! opcional: traza
        ! print *, iter, T_abs, T_sub, T_out, dx(1), dx(2), dx(3)
    end do

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

contains

    subroutine solve3x3(A, b, x)
        ! simple eliminacion gaussiana sin pivoteo parcial (stablea para casos normales)
        real(dp), intent(inout) :: A(3,3)
        real(dp), intent(in) :: b(3)
        real(dp), intent(out) :: x(3)
        real(dp) :: M(3,3), rhs(3)
        integer :: i,j,k
        real(dp) :: factor

        ! copiar para no modificar A externamente
        M = A
        rhs = b

        ! forward elimination
        do k = 1, 2
            if (abs(M(k,k)) < 1.0d-14) then
                print *, "Pivote casi cero en fila ", k, ". Abortando."
                stop
            end if
            do i = k+1, 3
                factor = M(i,k) / M(k,k)
                do j = k, 3
                    M(i,j) = M(i,j) - factor * M(k,j)
                end do
                rhs(i) = rhs(i) - factor * rhs(k)
            end do
        end do

        ! back substitution
        x(3) = rhs(3) / M(3,3)
        x(2) = (rhs(2) - M(2,3)*x(3)) / M(2,2)
        x(1) = (rhs(1) - M(1,2)*x(2) - M(1,3)*x(3)) / M(1,1)
    end subroutine solve3x3

end program Newton3Vars_HTF