program newton_ptc
    implicit none
    ! Declaración de variables
    real(8) :: Tptc, Tpv, Ta, Tsky
    real(8) :: Gghi, Aptc, alpha_ptc, eps_pv, eps_ptc
    real(8) :: sigma, Apv, Rconv
    real(8) :: f, df, Tnew
    real(8) :: tol
    integer :: iter, max_iter

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

    ! ==== Datos de ejemplo ====
    Gghi       = 800.0d0       ! W/m2
    Aptc       = 2.0d0         ! m2
    alpha_ptc  = 0.9d0
    eps_pv     = 0.85d0
    eps_ptc    = 0.9d0
    Apv        = 1.5d0         ! m2
    Tpv        = 330.0d0       ! K (~57°C)
    Ta         = 300.0d0       ! K (~27°C)
    Tsky       = 280.0d0       ! K (~7°C)
    Rconv      = 0.05d0        ! m2·K/W

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

    do iter = 1, max_iter
        ! Función f(Tptc)
        f = Gghi*Aptc*alpha_ptc + eps_pv*sigma*Apv*(Tpv**4 - Tptc**4) &
            - (Tptc - Ta)/Rconv - eps_ptc*sigma*Aptc*(Tptc**4 - Tsky**4)

        ! Derivada df/dTptc (analítica)
        df = -4.0d0*eps_pv*sigma*Apv*Tptc**3 - 1.0d0/Rconv &
             - 4.0d0*eps_ptc*sigma*Aptc*Tptc**3

        ! Actualización Newton–Raphson
        Tnew = Tptc - f/df

        ! Comprobar convergencia
        if (abs(Tnew - Tptc) < tol) then
            Tptc = Tnew
            exit
        endif

        Tptc = Tnew
    end do

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

end program newton_ptc