program ET_PTC
    
    
    ! Declaración de variables
    real(8) :: T_PTC, T_PVi, T_a, T_sky
    real(8) :: G_GHI, A_PTC, alpha_PTC, epsilon_PV, epsilon_PTC
    real(8) :: sigma, A_PV, R_conv_PTC
    real(8) :: f_T_PTC, df_T_PTC, T_PTCe
    real(8) :: tol
    integer :: iter, max_iter

    ! Constante de Stefan–Boltzmann (W/m2·K4)
    sigma        = 5.670374e-8

    ! ==== Datos de ejemplo ====
    G_GHI        = 1000.0d0      ! W/m2
    A_PTC        = 30.0d0         ! m2
    alpha_PTC    = 0.03d0
    epsilon_PV   = 0.2d0
    epsilon_PTC  = 0.3d0
    A_PV         = 1.2d0         ! m2
    T_PVi        = 340.5d0       ! K (~57°C)
    T_a          = 298.2d0       ! K (~27°C)
    T_sky        = 298.2d0       ! K (~7°C)
    R_conv_PTC   = 0.04336d0     ! m2·K/W

    ! ==== Configuración Newton–Raphson ====
    T_PTC = 320.0d0    ! Suposición inicial (~47°C)
    tol = 1.0d-6
    max_iter = 100



    do iter = 1, max_iter
        ! Función f(Tptc)
        f_T_PTC = G_GHI*A_PTC*alpha_PTC + epsilon_PV*sigma*A_PV*(T_PVi**4 - T_PTC**4) &
        - (T_PTC - T_a)/R_conv_PTC - epsilon_PTC*sigma*A_PTC*(T_PTC**4 - T_sky**4)

        ! Derivada df/dTptc (analítica)
        df_T_PTC = -4.0d0*epsilon_PV*sigma*A_PV*T_PTC**3 - 1.0d0/R_conv_PTC &
        - 4.0d0*epsilon_PTC*sigma*A_PTC*T_PTC**3

        ! Actualización Newton–Raphson
        T_PTCe = T_PTC - (f_T_PTC/df_T_PTC)

        ! Comprobar convergencia
        if (abs(T_PTCe - T_PTC) < tol) then
            T_PTC = T_PTCe
            exit
        endif

        T_PTC = T_PTCe
    end do



    ! Resultado
    print *, "Temperatura del PTC (K): ", T_PTC
    print *, "Temperatura del PTC (°C): ", T_PTC - 273.15d0
    print *, "Iteraciones: ", iter
    
end program ET_PTC