program combined_rad_conv
    implicit none

    ! Inputs
    double precision :: Ts, Tinf, Tsurr, eps, hc_dir, v, char_len, A
    integer :: hc_opt

    ! Temperatures in Kelvin
    double precision :: Ts_k, Tinf_k, Tsurr_k

    ! Calculated coefficients and heat rates
    double precision :: hc, hr, h_total
    double precision :: Q_rad, Q_conv, Q_total
    double precision :: pct_conv, pct_rad
    double precision :: delta_T_conv, delta_T_rad

    ! Constant
    double precision, parameter :: sigma = 5.670374d-8

    ! Read inputs from stdin
    read(*, *, iostat=hc_opt) Ts
    if (hc_opt /= 0) then
        print *, "ERROR: Failed to read Surface Temperature (Ts)"
        stop
    endif
    read(*, *) Tinf
    read(*, *) Tsurr
    read(*, *) eps
    read(*, *) hc_opt
    read(*, *) hc_dir
    read(*, *) v
    read(*, *) char_len
    read(*, *) A

    ! Validate Emissivity
    if (eps < 0.0d0 .or. eps > 1.0d0) then
        print *, "ERROR: Emissivity must be between 0.0 and 1.0"
        stop
    endif

    ! Validate Area
    if (A <= 0.0d0) then
        print *, "ERROR: Surface Area must be greater than 0"
        stop
    endif

    ! Convert temperatures to Kelvin
    Ts_k = Ts + 273.15d0
    Tinf_k = Tinf + 273.15d0
    Tsurr_k = Tsurr + 273.15d0

    ! Linearized radiation heat transfer coefficient hr
    hr = eps * sigma * (Ts_k + Tsurr_k) * (Ts_k**2 + Tsurr_k**2)
    Q_rad = hr * A * (Ts_k - Tsurr_k)

    ! Convective heat transfer coefficient hc
    select case (hc_opt)
    case (1) ! Direct input
        hc = hc_dir
    case (2) ! Natural convection, vertical plate
        delta_T_conv = abs(Ts - Tinf)
        hc = 1.42d0 * (delta_T_conv / max(0.001d0, char_len))**0.25d0
    case (3) ! Natural convection, horizontal cylinder
        delta_T_conv = abs(Ts - Tinf)
        hc = 1.32d0 * (delta_T_conv / max(0.001d0, char_len))**0.25d0
    case (4) ! Forced convection, air
        hc = 10.45d0 - v + 10.0d0 * sqrt(max(0.0d0, v))
    case default
        print *, "ERROR: Invalid convection option"
        stop
    end select

    ! Convective heat rate
    Q_conv = hc * A * (Ts_k - Tinf_k)

    ! Total heat rate
    Q_total = Q_conv + Q_rad

    ! Total heat transfer coefficient
    h_total = hc + hr

    ! Percentage calculation based on absolute heat rate contributions
    if (abs(Q_conv) + abs(Q_rad) > 1.0d-12) then
        pct_conv = (abs(Q_conv) / (abs(Q_conv) + abs(Q_rad))) * 100.0d0
        pct_rad = (abs(Q_rad) / (abs(Q_conv) + abs(Q_rad))) * 100.0d0
    else
        pct_conv = 50.0d0
        pct_rad = 50.0d0
    endif

    ! Print output in key-value format
    print *, "TS=", Ts
    print *, "TINF=", Tinf
    print *, "TSURR=", Tsurr
    print *, "EPS=", eps
    print *, "AREA=", A
    print *, "HC=", hc
    print *, "HR=", hr
    print *, "HTOTAL=", h_total
    print *, "Q_CONV=", Q_conv
    print *, "Q_RAD=", Q_rad
    print *, "Q_TOTAL=", Q_total
    print *, "PCT_CONV=", pct_conv
    print *, "PCT_RAD=", pct_rad

end program combined_rad_conv
