program Newton2Vars
    implicit none
    real(8) :: T_abs, T_sub, T_a, T_sky
    real(8) :: G_GHI, A_abs, alpha_abs, epsilon_abs
    real(8) :: sigma, R_conv_abs, R_cond_abs, R_cond_sub, R_conv_sub
    real(8) :: f1, f2, J11, J12, J21, J22
    real(8) :: dT_abs, dT_sub, detJ
    real(8) :: tol
    integer :: iter, max_iter

    sigma = 5.670374e-8

    ! ==== Datos de ejemplo ====
    G_GHI      = 1000
    A_abs      = 0.6
    alpha_abs  = 0.9
    epsilon_abs= 0.2
    T_a        = 298.2
    T_sky      = 298.2
    R_conv_abs = 0.0464
    R_cond_abs = 0.00002439
    R_cond_sub = 0.00004623
    R_conv_sub = 0.04

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

    do iter = 1, max_iter
        ! 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_abs - T_sub)/(R_cond_abs + R_cond_sub) &
            - (T_sub - T_a)/R_conv_sub

        ! Jacobiano
        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)
        J21 =  1.0d0/(R_cond_abs + R_cond_sub)
        J22 = -1.0d0/(R_cond_abs + R_cond_sub) - 1.0d0/R_conv_sub

        ! 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

        ! Comprobar 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