program combustion_aft
    implicit none
    integer :: fuel_type, iostat_val, i, n_sweep, iter
    double precision :: excess_air, T_inlet, P_kPa
    double precision :: C_at, H_at, O_at, M_fuel, LHV
    double precision :: a_stoich, air_stoich, AFR_stoich, AFR_actual
    double precision :: actual_O2, n_CO2, n_H2O, n_N2, n_O2_ex, n_total
    double precision :: y_CO2, y_H2O, y_N2, y_O2
    double precision :: T_ad, T_low, T_high, T_mid
    double precision :: H_react, H_prod, hf_fuel
    double precision :: cp_CO2, cp_H2O, cp_N2, cp_O2
    double precision :: hf_CO2, hf_H2O
    double precision :: P_H2O, T_dew
    double precision :: exc_sw, T_ad_sw, AFR_sw

    ! Average cp in J/(mol.K)
    cp_CO2 = 54.0d0
    cp_H2O = 41.0d0
    cp_N2  = 32.0d0
    cp_O2  = 34.0d0

    ! Enthalpies of formation in kJ/mol
    hf_CO2 = -393.5d0
    hf_H2O = -241.8d0

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) fuel_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid fuel type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) excess_air
    read(*,*,iostat=iostat_val) T_inlet
    read(*,*,iostat=iostat_val) P_kPa
    read(*,*,iostat=iostat_val) C_at
    read(*,*,iostat=iostat_val) H_at
    read(*,*,iostat=iostat_val) O_at
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (excess_air < 0.0d0) excess_air = 0.0d0
    if (T_inlet <= 0.0d0) T_inlet = 298.15d0
    if (P_kPa <= 0.0d0) P_kPa = 101.325d0

    ! ── Fuel properties ─────────────────────────────────────────
    if (fuel_type == 1) then
        C_at = 1.0d0; H_at = 4.0d0; O_at = 0.0d0
        M_fuel = 16.04d0; LHV = 50050.0d0
    else if (fuel_type == 2) then
        C_at = 3.0d0; H_at = 8.0d0; O_at = 0.0d0
        M_fuel = 44.10d0; LHV = 46357.0d0
    else if (fuel_type == 3) then
        C_at = 8.0d0; H_at = 18.0d0; O_at = 0.0d0
        M_fuel = 114.23d0; LHV = 44427.0d0
    else if (fuel_type == 4) then
        C_at = 0.0d0; H_at = 2.0d0; O_at = 0.0d0
        M_fuel = 2.016d0; LHV = 119960.0d0
    else
        ! Custom fuel
        if (C_at < 0.0d0) C_at = 1.0d0
        if (H_at < 0.0d0) H_at = 4.0d0
        if (O_at < 0.0d0) O_at = 0.0d0
        M_fuel = C_at*12.011d0 + H_at*1.008d0 + O_at*16.0d0
        if (M_fuel < 1.0d0) M_fuel = 16.04d0
        ! Estimate LHV from Dulong-like formula (kJ/kg)
        LHV = 33900.0d0*C_at*12.011d0/M_fuel + &
              142300.0d0*H_at*1.008d0/M_fuel - &
              9090.0d0*O_at*16.0d0/M_fuel
        if (LHV < 1000.0d0) LHV = 50000.0d0
    end if

    ! ── Stoichiometry ───────────────────────────────────────────
    ! CxHyOz + a(O2 + 3.76 N2) -> C CO2 + H/2 H2O + 3.76a N2
    ! a_stoich = C + H/4 - O_fuel/2
    a_stoich = C_at + H_at / 4.0d0 - O_at / 2.0d0
    if (a_stoich < 0.0d0) a_stoich = 0.001d0

    air_stoich = a_stoich / 0.21d0    ! moles air per mole fuel
    AFR_stoich = air_stoich * 28.97d0 / M_fuel   ! mass basis
    AFR_actual = AFR_stoich * (1.0d0 + excess_air / 100.0d0)

    ! Actual O2 supplied
    actual_O2 = a_stoich * (1.0d0 + excess_air / 100.0d0)

    ! Product moles per mole fuel
    n_CO2 = C_at
    n_H2O = H_at / 2.0d0
    n_N2  = 3.76d0 * actual_O2
    n_O2_ex = a_stoich * excess_air / 100.0d0
    n_total = n_CO2 + n_H2O + n_N2 + n_O2_ex

    ! Mole fractions
    y_CO2 = n_CO2 / max(n_total, 1.0d-10)
    y_H2O = n_H2O / max(n_total, 1.0d-10)
    y_N2  = n_N2  / max(n_total, 1.0d-10)
    y_O2  = n_O2_ex / max(n_total, 1.0d-10)

    ! ── Adiabatic flame temperature (bisection) ────────────────
    ! Enthalpy of formation of fuel from LHV:
    ! LHV = -hf_fuel + C*(-hf_CO2) + H/2*(-hf_H2O_g) [per mole fuel, roughly]
    ! Actually: hf_fuel = C*hf_CO2 + H/2*hf_H2O - (-LHV*M_fuel/1000)
    hf_fuel = n_CO2 * hf_CO2 + n_H2O * hf_H2O + LHV * M_fuel / 1000.0d0

    ! Energy balance: H_react = H_prod(T_ad)
    ! H_react = hf_fuel + actual_O2*(0) + n_N2_react*(0) + cp contributions at T_inlet
    ! H_react preheat = [actual_O2*cp_O2 + 3.76*actual_O2*cp_N2] * (T_inlet-298.15) / 1000
    H_react = hf_fuel + (actual_O2 * cp_O2 + 3.76d0*actual_O2 * cp_N2) * &
              (T_inlet - 298.15d0) / 1000.0d0

    ! H_prod(T) = n_CO2*hf_CO2 + n_H2O*hf_H2O + [n_CO2*cp_CO2 + n_H2O*cp_H2O + n_N2*cp_N2 + n_O2_ex*cp_O2]*(T-298.15)/1000
    ! Solve: H_react = H_prod(T_ad)
    ! => hf_fuel + preheat = n_CO2*hf_CO2 + n_H2O*hf_H2O + cp_prod_total*(T_ad-298.15)/1000
    ! => T_ad = 298.15 + 1000 * (hf_fuel + preheat - n_CO2*hf_CO2 - n_H2O*hf_H2O) / cp_prod_total

    T_low = 300.0d0
    T_high = 5000.0d0
    do iter = 1, 100
        T_mid = 0.5d0 * (T_low + T_high)
        H_prod = n_CO2*hf_CO2 + n_H2O*hf_H2O + &
                 (n_CO2*cp_CO2 + n_H2O*cp_H2O + n_N2*cp_N2 + n_O2_ex*cp_O2) * &
                 (T_mid - 298.15d0) / 1000.0d0
        if (H_prod < H_react) then
            T_low = T_mid
        else
            T_high = T_mid
        end if
        if (abs(T_high - T_low) < 0.1d0) exit
    end do
    T_ad = 0.5d0 * (T_low + T_high)

    ! ── Dew point ───────────────────────────────────────────────
    P_H2O = y_H2O * P_kPa   ! partial pressure of H2O in kPa
    ! Simplified Antoine for water: T_dew(K) ~ 273.15 + 100 * (P_H2O/101.325)^0.25
    ! More accurate simplified: T_dew = 273.15 + 45.0*log(P_H2O) + 1.0 (rough fit for low P)
    if (P_H2O > 0.0d0) then
        T_dew = 273.15d0 + 45.0d0 * log(max(P_H2O, 0.001d0)) + 7.0d0
    else
        T_dew = 273.15d0
    end if
    if (T_dew < 273.15d0) T_dew = 273.15d0

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   COMBUSTION & ADIABATIC FLAME TEMPERATURE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Fuel Type                 = ', fuel_type
    write(*,'(A,F10.2,A)')  '  Fuel Formula              : C', C_at, ''
    write(*,'(A,F10.2)')    '                              H', H_at
    if (O_at > 0.0d0) &
    write(*,'(A,F10.2)')    '                              O', O_at
    write(*,'(A,F10.2,A)')  '  Molar Mass M_fuel         = ', M_fuel, ' g/mol'
    write(*,'(A,F12.2,A)')  '  LHV                       = ', LHV, ' kJ/kg'
    write(*,'(A,F10.2,A)')  '  Excess Air                = ', excess_air, ' percent'
    write(*,'(A,F12.2,A)')  '  Inlet Temperature         = ', T_inlet, ' K'
    write(*,'(A,F12.4,A)')  '  Pressure                  = ', P_kPa, ' kPa'
    write(*,*)
    write(*,'(A)') '--- STOICHIOMETRY -------------------------------------------'
    write(*,'(A,F10.4)')    '  O2 stoichiometric (mol)   = ', a_stoich
    write(*,'(A,F10.4)')    '  Air stoich (mol/mol fuel) = ', air_stoich
    write(*,'(A,F10.4)')    '  AFR stoichiometric (mass) = ', AFR_stoich
    write(*,'(A,F10.4)')    '  AFR actual (mass)         = ', AFR_actual
    write(*,'(A,F10.4)')    '  Equivalence Ratio phi     = ', 1.0d0/(1.0d0+excess_air/100.0d0)
    write(*,*)
    write(*,'(A)') '--- FLUE GAS COMPOSITION (per mole fuel) -------------------'
    write(*,'(A,F10.4,A)')  '  CO2                       = ', n_CO2, ' mol'
    write(*,'(A,F10.4,A)')  '  H2O                       = ', n_H2O, ' mol'
    write(*,'(A,F10.4,A)')  '  N2                        = ', n_N2, ' mol'
    write(*,'(A,F10.4,A)')  '  O2 (excess)               = ', n_O2_ex, ' mol'
    write(*,'(A,F10.4,A)')  '  Total products            = ', n_total, ' mol'
    write(*,*)
    write(*,'(A)') '--- FLUE GAS MOLE FRACTIONS ---------------------------------'
    write(*,'(A,F10.4,A)')  '  y_CO2                     = ', y_CO2*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_H2O                     = ', y_H2O*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_N2                      = ', y_N2*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_O2                      = ', y_O2*100.0d0, ' percent'
    write(*,*)
    write(*,'(A)') '--- ADIABATIC FLAME TEMPERATURE -----------------------------'
    write(*,'(A,F12.2,A)')  '  T_ad                      = ', T_ad, ' K'
    write(*,'(A,F12.2,A)')  '  T_ad                      = ', T_ad - 273.15d0, ' C'
    write(*,*)
    write(*,'(A)') '--- DEW POINT -----------------------------------------------'
    write(*,'(A,F10.4,A)')  '  P_H2O (partial)           = ', P_H2O, ' kPa'
    write(*,'(A,F12.2,A)')  '  T_dew                     = ', T_dew, ' K'
    write(*,'(A,F12.2,A)')  '  T_dew                     = ', T_dew - 273.15d0, ' C'
    write(*,*)

    ! ── Sensitivity sweep: excess air 0% to 200% ──────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: T_AD VS EXCESS AIR -------------------------'
    write(*,'(A)') '  ExcessAir[%]  T_ad[K]       AFR_actual'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        exc_sw = dble(i-1) * 200.0d0 / dble(n_sweep - 1)
        ! Recompute for this excess air
        actual_O2 = a_stoich * (1.0d0 + exc_sw / 100.0d0)
        n_CO2 = C_at
        n_H2O = H_at / 2.0d0
        n_N2  = 3.76d0 * actual_O2
        n_O2_ex = a_stoich * exc_sw / 100.0d0
        n_total = n_CO2 + n_H2O + n_N2 + n_O2_ex

        H_react = hf_fuel + (actual_O2 * cp_O2 + 3.76d0*actual_O2 * cp_N2) * &
                  (T_inlet - 298.15d0) / 1000.0d0

        T_low = 300.0d0; T_high = 5000.0d0
        do iter = 1, 80
            T_mid = 0.5d0*(T_low+T_high)
            H_prod = n_CO2*hf_CO2 + n_H2O*hf_H2O + &
                     (n_CO2*cp_CO2+n_H2O*cp_H2O+n_N2*cp_N2+n_O2_ex*cp_O2) * &
                     (T_mid-298.15d0)/1000.0d0
            if (H_prod < H_react) then; T_low=T_mid; else; T_high=T_mid; end if
            if (abs(T_high-T_low)<0.1d0) exit
        end do
        T_ad_sw = 0.5d0*(T_low+T_high)
        AFR_sw = AFR_stoich * (1.0d0 + exc_sw/100.0d0)
        write(*,'(F10.2,4X,F12.2,4X,F10.4)') exc_sw, T_ad_sw, AFR_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Stoichiometric: CxHyOz + a(O2+3.76N2) -> CO2+H2O+N2'
    write(*,'(A)') '  a = C + H/4 - O_fuel/2'
    write(*,'(A)') '  AFR = a/0.21 * 28.97 / M_fuel (mass basis)'
    write(*,'(A)') '  T_ad by bisection on energy balance (constant cp)'
    write(*,'(A)') '  Dew point from partial pressure of H2O (simplified)'

end program combustion_aft