program solve_tem_conc
    implicit none
    real :: G_tot, A_conc, alpha_conc, sigma, epsilon_pv, A_pv, Tem_pv
    real :: h_conc, Tem_a, epsilon_conc, Tem_sky, Tem_conc
    real :: tol
    integer :: max_iter

    ! Asignar valores a las variables conocidas
    G_tot = 578
    A_conc = 94
    alpha_conc = 0.03
    sigma = 5.67e-8
    epsilon_pv = 0.2
    A_pv = 6
    Tem_pv = 355.45
    h_conc = 1.579
    Tem_a = 282
    epsilon_conc = 0.3
    Tem_sky = 282

    ! Parámetros para el método de Newton-Raphson
    tol = 1.0e-6
    max_iter = 100

    ! Estimación inicial de Tem_conc
    Tem_conc = 350.0

    ! Llamar a la función para resolver Tem_conc
    call newton_raphson(Tem_conc, tol, max_iter)

    ! Imprimir el resultado
    print *, 'La estimación de Tem_conc es: ', Tem_conc

contains

    function f(Tem_conc)
        real :: f, Tem_conc
        real :: q_solar_conc, q_rad_pv, q_conv_conc_front, q_conv_conc_back, q_rad_conc

        q_solar_conc = G_tot * A_conc * alpha_conc
        q_rad_pv = sigma * epsilon_pv * A_pv * (Tem_pv**4 - Tem_conc**4)
        q_conv_conc_front = A_conc * h_conc * (Tem_conc - Tem_a)
        q_conv_conc_back = A_conc * h_conc * (Tem_conc - Tem_a)
        q_rad_conc = sigma * epsilon_conc * A_conc * (Tem_conc**4 - Tem_sky**4)

        f = q_solar_conc + q_rad_pv - (q_conv_conc_front + q_conv_conc_back + q_rad_conc)
    end function f

    function df(Tem_conc)
        real :: df, Tem_conc
        real :: q_rad_pv, q_conv_conc, q_rad_conc

        q_rad_pv = sigma * epsilon_pv * A_pv * (-4.0 * Tem_conc**3)
        q_conv_conc = 2.0 * A_conc * h_conc
        q_rad_conc = sigma * epsilon_conc * A_conc * (4.0 * Tem_conc**3)

        df = q_rad_pv - (q_conv_conc + q_rad_conc)
    end function df

    subroutine newton_raphson(Tem_conc, tol, max_iter)
        real, intent(inout) :: Tem_conc
        real, intent(in) :: tol
        integer, intent(in) :: max_iter
        real :: f_val, df_val
        integer :: iter

        do iter = 1, max_iter
            f_val = f(Tem_conc)
            df_val = df(Tem_conc)

            if (abs(df_val) < tol) then
                print *, 'Error: La derivada es muy pequeña'
                stop
            endif

            Tem_conc = Tem_conc - f_val / df_val

            if (abs(f_val) < tol) then
                return
            endif
        end do

        print *, 'Advertencia: Número máximo de iteraciones alcanzado'
    end subroutine newton_raphson

end program solve_tem_conc