program fuel_cell_thermo
    implicit none
    integer :: cell_type, n_cells, iostat_val, i, n_sweep
    double precision :: T_cell, P_cell, i_dens, A_active
    double precision :: P_H2, P_O2, alpha_a, i0_a, r_ohm, i_lim
    double precision :: n_e, dG0, dH0, dS0, E0_rev, M_fuel
    double precision :: F_const, R_gas
    double precision :: E_rev_T, E_Nernst, V_act, V_ohm, V_conc, V_cell
    double precision :: P_density, P_cell_W, P_stack
    double precision :: eta_thermo, eta_voltage, eta_overall
    double precision :: Q_heat, mdot_fuel
    double precision :: i_sw, V_sw, P_sw
    character(len=50) :: cell_name, fuel_name

    F_const = 96485.0d0   ! C/mol
    R_gas = 8.314462d0     ! J/(mol.K)

    read(*,*,iostat=iostat_val) cell_type
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid cell type.'; stop; end if
    read(*,*,iostat=iostat_val) T_cell
    read(*,*,iostat=iostat_val) P_cell
    read(*,*,iostat=iostat_val) i_dens
    read(*,*,iostat=iostat_val) A_active
    read(*,*,iostat=iostat_val) n_cells
    read(*,*,iostat=iostat_val) P_H2
    read(*,*,iostat=iostat_val) P_O2
    read(*,*,iostat=iostat_val) alpha_a
    read(*,*,iostat=iostat_val) i0_a
    read(*,*,iostat=iostat_val) r_ohm
    read(*,*,iostat=iostat_val) i_lim
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if

    if (T_cell <= 0.0d0) T_cell = 353.15d0
    if (P_cell <= 0.0d0) P_cell = 1.0d0
    if (i_dens <= 0.0d0) i_dens = 0.5d0
    if (A_active <= 0.0d0) A_active = 100.0d0
    if (n_cells < 1) n_cells = 1
    if (P_H2 <= 0.0d0) P_H2 = 1.0d0
    if (P_O2 <= 0.0d0) P_O2 = 0.21d0
    if (alpha_a <= 0.0d0) alpha_a = 0.5d0
    if (i0_a <= 0.0d0) i0_a = 1.0d-4
    if (r_ohm <= 0.0d0) r_ohm = 0.15d0
    if (i_lim <= 0.0d0) i_lim = 2.0d0

    ! ── Cell reaction data ─────────────────────────────────────
    select case(cell_type)
    case(1)
        cell_name = 'PEM Fuel Cell (H2/O2)'
        fuel_name = 'Hydrogen H2'
        n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
        M_fuel = 2.016d0
    case(2)
        cell_name = 'Solid Oxide Fuel Cell (SOFC, H2/O2)'
        fuel_name = 'Hydrogen H2'
        n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
        M_fuel = 2.016d0
        if (T_cell < 873.15d0) T_cell = 1073.15d0  ! SOFC operates at high T
    case(3)
        cell_name = 'Direct Methanol FC (DMFC)'
        fuel_name = 'Methanol CH3OH'
        n_e = 6.0d0; dG0 = -702.5d0; dH0 = -726.6d0; E0_rev = 1.213d0
        M_fuel = 32.04d0
    case default
        cell_name = 'Alkaline Fuel Cell (AFC, H2/O2)'
        fuel_name = 'Hydrogen H2'
        n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
        M_fuel = 2.016d0
        cell_type = 4
    end select

    ! ── Thermodynamic calculations ─────────────────────────────
    ! Entropy change
    dS0 = (dH0 - dG0) * 1000.0d0 / 298.15d0   ! J/(mol.K)

    ! Reversible voltage at temperature T
    E_rev_T = E0_rev + dS0 / (n_e * F_const) * (T_cell - 298.15d0)

    ! Nernst equation (partial pressure correction)
    ! For H2 + 0.5O2 -> H2O:
    ! E_Nernst = E_rev(T) + (RT/(n_e*F))*ln(P_H2 * P_O2^0.5)
    E_Nernst = E_rev_T + (R_gas * T_cell / (n_e * F_const)) * &
               log(P_H2 * sqrt(P_O2))

    ! ── Polarization losses ────────────────────────────────────
    ! Activation overpotential (Butler-Volmer / Tafel)
    if (i_dens > 0.0d0 .and. i0_a > 0.0d0) then
        V_act = (R_gas * T_cell / (alpha_a * n_e * F_const)) * &
                log(i_dens / i0_a + sqrt((i_dens/i0_a)**2 + 1.0d0))
    else
        V_act = 0.0d0
    end if

    ! Ohmic losses
    V_ohm = i_dens * r_ohm

    ! Concentration (mass transport) losses
    if (i_dens < i_lim .and. i_dens > 0.0d0) then
        V_conc = -(R_gas * T_cell / (n_e * F_const)) * log(1.0d0 - i_dens / i_lim)
    else
        V_conc = 0.5d0   ! saturated
    end if

    ! Cell voltage
    V_cell = E_Nernst - V_act - V_ohm - V_conc
    if (V_cell < 0.0d0) V_cell = 0.0d0

    ! Power density
    P_density = V_cell * i_dens   ! W/cm2

    ! Cell and stack power
    P_cell_W = V_cell * i_dens * A_active   ! W
    P_stack = P_cell_W * dble(n_cells)       ! W

    ! Efficiencies
    eta_thermo = abs(dG0 / dH0)   ! max thermodynamic efficiency
    if (E_Nernst > 1.0d-10) then
        eta_voltage = V_cell / E_Nernst
    else
        eta_voltage = 0.0d0
    end if
    eta_overall = eta_voltage * eta_thermo

    ! Heat generation
    Q_heat = (E_Nernst - V_cell) * i_dens * A_active * dble(n_cells)   ! W

    ! Fuel consumption
    mdot_fuel = i_dens * A_active * dble(n_cells) * M_fuel / (n_e * F_const)   ! g/s

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   FUEL CELL THERMODYNAMICS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Cell Type                 = ', trim(cell_name)
    write(*,'(A,A)')        '  Fuel                      = ', trim(fuel_name)
    write(*,'(A,F12.2,A)')  '  Cell Temperature          = ', T_cell, ' K'
    write(*,'(A,F12.2,A)')  '  Cell Temperature          = ', T_cell-273.15d0, ' C'
    write(*,'(A,F10.4,A)')  '  Current Density           = ', i_dens, ' A/cm2'
    write(*,'(A,F12.2,A)')  '  Active Area               = ', A_active, ' cm2'
    write(*,'(A,I8)')       '  Number of Cells           = ', n_cells
    write(*,'(A,F10.4,A)')  '  P_H2                      = ', P_H2, ' atm'
    write(*,'(A,F10.4,A)')  '  P_O2                      = ', P_O2, ' atm'
    write(*,*)
    write(*,'(A)') '--- ELECTROCHEMICAL PARAMETERS ------------------------------'
    write(*,'(A,F10.4)')    '  Transfer Coefficient a    = ', alpha_a
    write(*,'(A,ES12.4,A)') '  Exchange Current i0       = ', i0_a, ' A/cm2'
    write(*,'(A,F10.4,A)')  '  Ohmic Resistance r        = ', r_ohm, ' ohm.cm2'
    write(*,'(A,F10.4,A)')  '  Limiting Current i_lim    = ', i_lim, ' A/cm2'
    write(*,*)
    write(*,'(A)') '--- THERMODYNAMICS ------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Delta_G0                  = ', dG0, ' kJ/mol'
    write(*,'(A,F12.2,A)')  '  Delta_H0                  = ', dH0, ' kJ/mol'
    write(*,'(A,F12.4,A)')  '  Delta_S0                  = ', dS0, ' J/(mol.K)'
    write(*,'(A,F10.6)')    '  n_electrons               = ', n_e
    write(*,'(A,F10.6,A)')  '  E0_rev (25C)              = ', E0_rev, ' V'
    write(*,'(A,F10.6,A)')  '  E_rev(T)                  = ', E_rev_T, ' V'
    write(*,'(A,F10.6,A)')  '  E_Nernst                  = ', E_Nernst, ' V'
    write(*,*)
    write(*,'(A)') '--- POLARIZATION LOSSES -------------------------------------'
    write(*,'(A,F10.6,A)')  '  Activation V_act          = ', V_act, ' V'
    write(*,'(A,F10.6,A)')  '  Ohmic V_ohm               = ', V_ohm, ' V'
    write(*,'(A,F10.6,A)')  '  Concentration V_conc      = ', V_conc, ' V'
    write(*,'(A,F10.6,A)')  '  Total Losses              = ', V_act+V_ohm+V_conc, ' V'
    write(*,*)
    write(*,'(A)') '--- CELL PERFORMANCE ----------------------------------------'
    write(*,'(A,F10.6,A)')  '  Cell Voltage V_cell       = ', V_cell, ' V'
    write(*,'(A,F10.6,A)')  '  Power Density             = ', P_density, ' W/cm2'
    write(*,'(A,F12.4,A)')  '  Cell Power                = ', P_cell_W, ' W'
    write(*,'(A,F12.4,A)')  '  Stack Power               = ', P_stack, ' W'
    write(*,'(A,F12.2,A)')  '  Stack Power               = ', P_stack/1000.0d0, ' kW'
    write(*,*)
    write(*,'(A)') '--- EFFICIENCY ----------------------------------------------'
    write(*,'(A,F10.4)')    '  Thermo Efficiency (max)   = ', eta_thermo
    write(*,'(A,F10.2,A)')  '  Thermo Efficiency         = ', eta_thermo*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  Voltage Efficiency        = ', eta_voltage
    write(*,'(A,F10.2,A)')  '  Voltage Efficiency        = ', eta_voltage*100.0d0, ' percent'
    write(*,'(A,F10.4)')    '  Overall Efficiency        = ', eta_overall
    write(*,'(A,F10.2,A)')  '  Overall Efficiency        = ', eta_overall*100.0d0, ' percent'
    write(*,*)
    write(*,'(A)') '--- HEAT & FUEL CONSUMPTION ---------------------------------'
    write(*,'(A,F12.4,A)')  '  Heat Generation           = ', Q_heat, ' W'
    write(*,'(A,F12.6,A)')  '  Fuel Consumption          = ', mdot_fuel, ' g/s'
    write(*,'(A,F12.4,A)')  '  Fuel Consumption          = ', mdot_fuel*3600.0d0, ' g/h'
    write(*,*)

    ! ── Polarization curve sweep ──────────────────────────────
    n_sweep = 50
    write(*,'(A)') '--- POLARIZATION CURVE --------------------------------------'
    write(*,'(A)') '  i[A/cm2]      V_cell[V]     P_dens[W/cm2]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        i_sw = 0.01d0 + dble(i-1) * (i_lim*0.95d0 - 0.01d0) / dble(n_sweep-1)
        if (i_sw > 0.0d0 .and. i0_a > 0.0d0) then
            V_act = (R_gas*T_cell/(alpha_a*n_e*F_const)) * &
                    log(i_sw/i0_a + sqrt((i_sw/i0_a)**2+1.0d0))
        else
            V_act = 0.0d0
        end if
        V_ohm = i_sw * r_ohm
        if (i_sw < i_lim) then
            V_conc = -(R_gas*T_cell/(n_e*F_const))*log(1.0d0-i_sw/i_lim)
        else
            V_conc = 0.5d0
        end if
        V_sw = E_Nernst - V_act - V_ohm - V_conc
        if (V_sw < 0.0d0) V_sw = 0.0d0
        P_sw = V_sw * i_sw
        write(*,'(F10.4,4X,F10.6,4X,F10.6)') i_sw, V_sw, P_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Nernst: E = E0 + (dS/(nF))(T-298) + (RT/(nF))ln(PH2*PO2^0.5)'
    write(*,'(A)') '  Activation: V_act = (RT/(a*n*F))*arcsinh(i/(2*i0))'
    write(*,'(A)') '  Ohmic: V_ohm = i * R_area'
    write(*,'(A)') '  Concentration: V_conc = -(RT/(nF))*ln(1 - i/i_lim)'
    write(*,'(A)') '  Efficiency: eta = (dG/dH) * (V_cell/E_Nernst)'

end program fuel_cell_thermo
