program real_gas_eos
    implicit none
    integer :: eos_type, i, iostat_val
    double precision :: T, P, Tc, Pc, omega, Rgas, Mmol
    double precision :: a_eos, b_eos, kappa, alpha_pr
    double precision :: v, Z, v_ideal, Z_ideal
    double precision :: dh_dep, ds_dep, dg_dep, fugacity, phi_fug
    double precision :: P_i, v_i, Z_i, T_i
    double precision, parameter :: Ru = 8.314462d0   ! J/(mol K)
    character(len=40) :: eos_name

    read(*,*,iostat=iostat_val) eos_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid EOS type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T
    read(*,*,iostat=iostat_val) P
    read(*,*,iostat=iostat_val) Tc
    read(*,*,iostat=iostat_val) Pc
    read(*,*,iostat=iostat_val) omega
    read(*,*,iostat=iostat_val) Mmol
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all EOS inputs.'
        stop
    end if
    if (T <= 0.0d0 .or. P <= 0.0d0 .or. Tc <= 0.0d0 .or. Pc <= 0.0d0) then
        write(*,*) 'ERROR: T, P, Tc, Pc must be positive.'
        stop
    end if
    if (Mmol <= 0.0d0) Mmol = 28.97d0
    Rgas = Ru / (Mmol * 1.0d-3)   ! J/(kg K) if Mmol in g/mol

    ! Ideal gas
    v_ideal = Ru * T / P           ! m3/mol
    Z_ideal = 1.0d0

    select case(eos_type)
    case(1)
        eos_name = 'van der Waals'
        a_eos = 27.0d0 * Ru**2 * Tc**2 / (64.0d0 * Pc)
        b_eos = Ru * Tc / (8.0d0 * Pc)
        call solve_cubic_vdw(T, P, a_eos, b_eos, Ru, v)
    case(2)
        eos_name = 'Redlich-Kwong'
        a_eos = 0.42748d0 * Ru**2 * Tc**2.5d0 / Pc
        b_eos = 0.08664d0 * Ru * Tc / Pc
        call solve_cubic_rk(T, P, a_eos, b_eos, Ru, v)
    case(3)
        eos_name = 'Peng-Robinson'
        kappa = 0.37464d0 + 1.54226d0*omega - 0.26992d0*omega**2
        alpha_pr = (1.0d0 + kappa*(1.0d0 - sqrt(T/Tc)))**2
        a_eos = 0.45724d0 * Ru**2 * Tc**2 / Pc * alpha_pr
        b_eos = 0.07780d0 * Ru * Tc / Pc
        call solve_cubic_pr(T, P, a_eos, b_eos, Ru, v)
    case default
        write(*,*) 'ERROR: EOS type must be 1 vdW, 2 RK, or 3 PR.'
        stop
    end select

    Z = P * v / (Ru * T)

    ! Departure functions (molar basis, J/mol)
    call departure_functions(eos_type, T, P, v, a_eos, b_eos, Ru, Tc, &
                             dh_dep, ds_dep, dg_dep, fugacity, phi_fug)

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   REAL GAS — EQUATIONS OF STATE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Equation of State         = ', trim(eos_name)
    write(*,'(A,ES12.4,A)') '  Temperature T             = ', T, ' K'
    write(*,'(A,ES12.4,A)') '  Pressure P                = ', P, ' Pa'
    write(*,'(A,ES12.4,A)') '  Critical Temperature Tc   = ', Tc, ' K'
    write(*,'(A,ES12.4,A)') '  Critical Pressure Pc      = ', Pc, ' Pa'
    write(*,'(A,ES12.4)')   '  Acentric Factor omega     = ', omega
    write(*,'(A,ES12.4,A)') '  Molar Mass                = ', Mmol, ' g/mol'
    write(*,'(A,ES12.4)')   '  Reduced Temperature Tr    = ', T/Tc
    write(*,'(A,ES12.4)')   '  Reduced Pressure Pr       = ', P/Pc
    write(*,*)
    write(*,'(A)') '--- EOS PARAMETERS ------------------------------------------'
    write(*,'(A,ES12.4)')   '  a (attraction)            = ', a_eos
    write(*,'(A,ES12.4)')   '  b (co-volume)             = ', b_eos
    write(*,*)
    write(*,'(A)') '--- STATE RESULTS -------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Molar Volume v            = ', v, ' m3/mol'
    write(*,'(A,ES12.4,A)') '  Ideal Gas Volume          = ', v_ideal, ' m3/mol'
    write(*,'(A,ES12.6)')   '  Compressibility Factor Z  = ', Z
    write(*,'(A,ES12.4,A)') '  Specific Volume           = ', v/(Mmol*1.0d-3), ' m3/kg'
    write(*,'(A,ES12.4,A)') '  Density                   = ', (Mmol*1.0d-3)/v, ' kg/m3'
    write(*,*)
    write(*,'(A)') '--- DEPARTURE FUNCTIONS (molar) -----------------------------'
    write(*,'(A,ES12.4,A)') '  Enthalpy Departure dh     = ', dh_dep, ' J/mol'
    write(*,'(A,ES12.4,A)') '  Entropy Departure ds      = ', ds_dep, ' J/(mol.K)'
    write(*,'(A,ES12.4,A)') '  Gibbs Departure dg        = ', dg_dep, ' J/mol'
    write(*,'(A,ES12.4,A)') '  Fugacity                  = ', fugacity, ' Pa'
    write(*,'(A,ES12.6)')   '  Fugacity Coefficient phi  = ', phi_fug
    write(*,*)

    ! Z vs P sweep (isothermal)
    write(*,'(A)') '--- Z VS PRESSURE SWEEP (isothermal) ------------------------'
    write(*,'(A)') '  P[Pa]         Z             v[m3/mol]     rho[kg/m3]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 60
        P_i = P * 0.01d0 * (200.0d0)**(dble(i-1)/59.0d0)
        select case(eos_type)
        case(1); call solve_cubic_vdw(T, P_i, a_eos, b_eos, Ru, v_i)
        case(2); call solve_cubic_rk(T, P_i, a_eos, b_eos, Ru, v_i)
        case(3); call solve_cubic_pr(T, P_i, a_eos, b_eos, Ru, v_i)
        end select
        Z_i = P_i * v_i / (Ru * T)
        write(*,'(ES12.4,2X,F12.6,2X,ES12.4,2X,ES12.4)') &
            P_i, Z_i, v_i, (Mmol*1.0d-3)/v_i
    end do
    write(*,*)

    ! Z vs T sweep (isobaric)
    write(*,'(A)') '--- Z VS TEMPERATURE SWEEP (isobaric) -----------------------'
    write(*,'(A)') '  T[K]          Z             v[m3/mol]     Tr'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 50
        T_i = Tc * 0.5d0 + (Tc*3.0d0 - Tc*0.5d0)*dble(i-1)/49.0d0
        ! Recompute a for PR (temperature-dependent alpha)
        if (eos_type == 3) then
            alpha_pr = (1.0d0 + kappa*(1.0d0 - sqrt(T_i/Tc)))**2
            a_eos = 0.45724d0 * Ru**2 * Tc**2 / Pc * alpha_pr
        end if
        if (eos_type == 2) then
            a_eos = 0.42748d0 * Ru**2 * Tc**2.5d0 / Pc
        end if
        select case(eos_type)
        case(1); call solve_cubic_vdw(T_i, P, a_eos, b_eos, Ru, v_i)
        case(2); call solve_cubic_rk(T_i, P, a_eos, b_eos, Ru, v_i)
        case(3); call solve_cubic_pr(T_i, P, a_eos, b_eos, Ru, v_i)
        end select
        Z_i = P * v_i / (Ru * T_i)
        write(*,'(F12.2,2X,F12.6,2X,ES12.4,2X,F10.4)') T_i, Z_i, v_i, T_i/Tc
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  vdW: P = RT/(v-b) - a/v^2.'
    write(*,'(A)') '  RK:  P = RT/(v-b) - a/(T^0.5 v(v+b)).'
    write(*,'(A)') '  PR:  P = RT/(v-b) - a alpha/(v(v+b)+b(v-b)).'
    write(*,'(A)') '  Z = Pv/(RT).'

contains

    subroutine solve_cubic_vdw(Tin, Pin, ain, bin, R, vout)
        implicit none
        double precision, intent(in) :: Tin, Pin, ain, bin, R
        double precision, intent(out) :: vout
        double precision :: v_guess, Pv, dPdv, dv
        integer :: it
        v_guess = R*Tin/Pin
        do it = 1, 200
            Pv = R*Tin/(v_guess-bin) - ain/v_guess**2
            dPdv = -R*Tin/(v_guess-bin)**2 + 2.0d0*ain/v_guess**3
            dv = -(Pv - Pin)/dPdv
            v_guess = v_guess + dv
            if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
            if (abs(dv) < 1.0d-15*abs(v_guess)) exit
        end do
        vout = v_guess
    end subroutine solve_cubic_vdw

    subroutine solve_cubic_rk(Tin, Pin, ain, bin, R, vout)
        implicit none
        double precision, intent(in) :: Tin, Pin, ain, bin, R
        double precision, intent(out) :: vout
        double precision :: v_guess, Pv, dPdv, dv, sqT
        integer :: it
        sqT = sqrt(Tin)
        v_guess = R*Tin/Pin
        do it = 1, 200
            Pv = R*Tin/(v_guess-bin) - ain/(sqT*v_guess*(v_guess+bin))
            dPdv = -R*Tin/(v_guess-bin)**2 + ain/(sqT)* &
                   (2.0d0*v_guess+bin)/(v_guess*(v_guess+bin))**2
            dv = -(Pv - Pin)/dPdv
            v_guess = v_guess + dv
            if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
            if (abs(dv) < 1.0d-15*abs(v_guess)) exit
        end do
        vout = v_guess
    end subroutine solve_cubic_rk

    subroutine solve_cubic_pr(Tin, Pin, ain, bin, R, vout)
        implicit none
        double precision, intent(in) :: Tin, Pin, ain, bin, R
        double precision, intent(out) :: vout
        double precision :: v_guess, Pv, dPdv, dv, denom
        integer :: it
        v_guess = R*Tin/Pin
        do it = 1, 200
            denom = v_guess*(v_guess+bin)+bin*(v_guess-bin)
            Pv = R*Tin/(v_guess-bin) - ain/denom
            dPdv = -R*Tin/(v_guess-bin)**2 + ain*(2.0d0*v_guess+2.0d0*bin-bin)/ &
                   denom**2 * (2.0d0*v_guess+2.0d0*bin)
            ! Numerical derivative fallback
            dPdv = -R*Tin/(v_guess-bin)**2 + &
                   ain*(2.0d0*v_guess)/(denom**2)
            dv = -(Pv - Pin)/dPdv
            if (abs(dv) > 0.5d0*v_guess) dv = sign(0.5d0*v_guess, dv)
            v_guess = v_guess + dv
            if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
            if (abs(dv) < 1.0d-15*abs(v_guess)) exit
        end do
        vout = v_guess
    end subroutine solve_cubic_pr

    subroutine departure_functions(etype, Tin, Pin, vin, ain, bin, R, Tcin, &
                                   dh, ds, dg, fug, phi)
        implicit none
        integer, intent(in) :: etype
        double precision, intent(in) :: Tin, Pin, vin, ain, bin, R, Tcin
        double precision, intent(out) :: dh, ds, dg, fug, phi
        double precision :: Zval, lnphi

        Zval = Pin*vin/(R*Tin)

        select case(etype)
        case(1)
            ! van der Waals departure
            dh = R*Tin*(Zval - 1.0d0) - 2.0d0*ain/vin + Pin*vin - R*Tin
            dh = Pin*vin - R*Tin - ain/vin
            ds = R*log(max((vin-bin)*Pin/(R*Tin), 1.0d-30))
            lnphi = bin/(vin-bin) - 2.0d0*ain/(R*Tin*vin) - log(Zval)
        case(2)
            ! Redlich-Kwong departure
            dh = R*Tin*(Zval-1.0d0) - 1.5d0*ain/(bin*sqrt(Tin))* &
                 log((vin+bin)/vin)
            ds = R*log(max(Zval*(1.0d0-bin/vin),1.0d-30)) - &
                 0.5d0*ain/(bin*Tin*sqrt(Tin))*log((vin+bin)/vin)
            lnphi = Zval - 1.0d0 - log(max(Zval*(1.0d0-bin/vin),1.0d-30)) - &
                     ain/(bin*R*sqrt(Tin)*Tin)*log((vin+bin)/vin)
        case(3)
            ! Peng-Robinson departure
            dh = R*Tin*(Zval-1.0d0) + (Tin*0.0d0 - ain)/(2.0d0*sqrt(2.0d0)*bin)* &
                 log((vin+(1.0d0+sqrt(2.0d0))*bin)/(vin+(1.0d0-sqrt(2.0d0))*bin))
            ds = R*log(max(Zval*(vin-bin)/vin, 1.0d-30))
            lnphi = Zval - 1.0d0 - log(max(Zval-bin*Pin/(R*Tin),1.0d-30)) - &
                     ain/(2.0d0*sqrt(2.0d0)*bin*R*Tin)* &
                     log((vin+(1.0d0+sqrt(2.0d0))*bin)/ &
                         (vin+(1.0d0-sqrt(2.0d0))*bin))
        end select

        phi = exp(lnphi)
        fug = phi * Pin
        dg = R*Tin*lnphi
    end subroutine departure_functions

end program real_gas_eos