program isentropic_polytropic
    implicit none
    integer :: i, iostat_val, n_sweep
    double precision :: P1, P2, T1, V1, poly_n, R_gas, gamma_val, mass
    double precision :: V2, T2, W, Q, dU, dH, dS
    double precision :: cv, cp
    double precision :: n_sw, V2_sw, T2_sw, W_sw, Q_sw
    character(len=40) :: process_name

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) P1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid P1 input.'
        stop
    end if
    read(*,*,iostat=iostat_val) P2
    read(*,*,iostat=iostat_val) T1
    read(*,*,iostat=iostat_val) V1
    read(*,*,iostat=iostat_val) poly_n
    read(*,*,iostat=iostat_val) R_gas
    read(*,*,iostat=iostat_val) gamma_val
    read(*,*,iostat=iostat_val) mass
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (P1 <= 0.0d0) then
        write(*,*) 'ERROR: P1 must be positive.'
        stop
    end if
    if (P2 <= 0.0d0) then
        write(*,*) 'ERROR: P2 must be positive.'
        stop
    end if
    if (T1 <= 0.0d0) T1 = 300.0d0
    if (V1 <= 0.0d0) V1 = 0.1d0
    if (mass <= 0.0d0) mass = 1.0d0
    if (gamma_val <= 1.0d0) gamma_val = 1.4d0
    if (R_gas <= 0.0d0) R_gas = 287.0d0

    ! ── Gas properties ──────────────────────────────────────────
    cv = R_gas / (gamma_val - 1.0d0)       ! J/(kg·K)
    cp = gamma_val * cv                      ! J/(kg·K)

    ! ── Process identification ──────────────────────────────────
    if (abs(poly_n) < 1.0d-6) then
        process_name = 'Isobaric (n = 0)'
    else if (abs(poly_n - 1.0d0) < 1.0d-6) then
        process_name = 'Isothermal (n = 1)'
    else if (abs(poly_n - gamma_val) < 0.02d0) then
        process_name = 'Isentropic (n = gamma)'
    else if (poly_n > 20.0d0) then
        process_name = 'Isochoric (n -> infinity)'
    else
        process_name = 'Polytropic (general)'
    end if

    ! ── Core calculations ───────────────────────────────────────
    ! Handle special cases for n
    if (poly_n > 20.0d0) then
        ! Isochoric: V = const
        V2 = V1
        T2 = T1 * (P2 / P1)
        W = 0.0d0
    else if (abs(poly_n) < 1.0d-6) then
        ! Isobaric: P = const => P2 = P1 (force)
        V2 = V1 * (T1 * P2) / (T1 * P1)  ! from PV=mRT
        ! Actually for isobaric: V2/V1 = T2/T1, need another relation
        ! With n=0: PV^0 = C => P = C, so P2 = P1
        ! T2 = T1 * V2/V1, but V2 = V1*(P1/P2)^(1/0) is undefined
        ! Use ideal gas: T2/T1 = (P2/P1)^((n-1)/n) with n->0: T2/T1 = (P2/P1)^(-inf)
        ! For n=0, P=const so P2=P1. Use V2 from ideal gas
        T2 = T1 * (P2 / P1)  ! simplified
        V2 = V1 * T2 / T1
        W = P1 * (V2 - V1)   ! kPa * m^3 = kJ
    else if (abs(poly_n - 1.0d0) < 1.0d-6) then
        ! Isothermal: T = const
        T2 = T1
        V2 = V1 * (P1 / P2)
        W = P1 * V1 * log(V2 / V1)   ! kPa * m^3 = kJ
    else
        ! General polytropic
        V2 = V1 * (P1 / P2) ** (1.0d0 / poly_n)
        T2 = T1 * (P2 / P1) ** ((poly_n - 1.0d0) / poly_n)
        W = (P2 * V2 - P1 * V1) / (1.0d0 - poly_n)   ! kPa*m^3 = kJ
    end if

    ! Internal energy change (kJ)
    dU = mass * cv * (T2 - T1) / 1000.0d0

    ! Enthalpy change (kJ)
    dH = mass * cp * (T2 - T1) / 1000.0d0

    ! Heat transfer from first law: Q = dU + W (all in kJ)
    Q = dU + W

    ! Entropy change for ideal gas (kJ/K)
    if (T2 > 0.0d0 .and. T1 > 0.0d0 .and. V2 > 0.0d0 .and. V1 > 0.0d0) then
        dS = mass * (cv * log(T2/T1) + R_gas * log(V2/V1)) / 1000.0d0
    else
        dS = 0.0d0
    end if

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ISENTROPIC & POLYTROPIC PROCESS ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,F12.4,A)')  '  Initial Pressure P1       = ', P1, ' kPa'
    write(*,'(A,F12.4,A)')  '  Final Pressure P2         = ', P2, ' kPa'
    write(*,'(A,F12.2,A)')  '  Initial Temperature T1    = ', T1, ' K'
    write(*,'(A,ES12.4,A)') '  Initial Volume V1         = ', V1, ' m^3'
    write(*,'(A,F12.4)')    '  Polytropic Index n        = ', poly_n
    write(*,'(A,F12.2,A)')  '  Specific Gas Constant R   = ', R_gas, ' J/(kg.K)'
    write(*,'(A,F12.4)')    '  Gamma (cp/cv)             = ', gamma_val
    write(*,'(A,F12.4,A)')  '  Mass                      = ', mass, ' kg'
    write(*,*)
    write(*,'(A)') '--- GAS PROPERTIES ------------------------------------------'
    write(*,'(A,F12.4,A)')  '  cv                        = ', cv, ' J/(kg.K)'
    write(*,'(A,F12.4,A)')  '  cp                        = ', cp, ' J/(kg.K)'
    write(*,*)
    write(*,'(A)') '--- PROCESS IDENTIFICATION ----------------------------------'
    write(*,'(A,A)')        '  Process Type              = ', trim(process_name)
    write(*,*)
    write(*,'(A)') '--- FINAL STATE ---------------------------------------------'
    write(*,'(A,F12.2,A)')  '  Final Temperature T2      = ', T2, ' K'
    write(*,'(A,ES12.4,A)') '  Final Volume V2           = ', V2, ' m^3'
    write(*,'(A,F12.4,A)')  '  Final Pressure P2         = ', P2, ' kPa'
    write(*,*)
    write(*,'(A)') '--- ENERGY & ENTROPY ----------------------------------------'
    write(*,'(A,F12.4,A)')  '  Work W                    = ', W, ' kJ'
    write(*,'(A,F12.4,A)')  '  Heat Transfer Q           = ', Q, ' kJ'
    write(*,'(A,F12.4,A)')  '  Internal Energy Change dU = ', dU, ' kJ'
    write(*,'(A,F12.4,A)')  '  Enthalpy Change dH        = ', dH, ' kJ'
    write(*,'(A,ES14.6,A)') '  Entropy Change dS         = ', dS, ' kJ/K'
    write(*,*)
    write(*,'(A)') '--- SPECIFIC VALUES (per kg) --------------------------------'
    if (mass > 1.0d-10) then
    write(*,'(A,F12.4,A)')  '  w (specific work)         = ', W/mass, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  q (specific heat)         = ', Q/mass, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  du                        = ', dU/mass, ' kJ/kg'
    write(*,'(A,F12.4,A)')  '  dh                        = ', dH/mass, ' kJ/kg'
    write(*,'(A,ES14.6,A)') '  ds                        = ', dS/mass, ' kJ/(kg.K)'
    end if
    write(*,*)

    ! ── Sensitivity sweep: n from 0.2 to 3.0 ───────────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: POLYTROPIC INDEX SWEEP ---------------------'
    write(*,'(A)') '  n           T2[K]         W[kJ]         Q[kJ]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        n_sw = 0.2d0 + dble(i-1) * (3.0d0 - 0.2d0) / dble(n_sweep - 1)
        if (abs(n_sw - 1.0d0) < 0.02d0) then
            T2_sw = T1
            V2_sw = V1 * (P1 / P2)
            W_sw = P1 * V1 * log(V2_sw / V1)
        else
            V2_sw = V1 * (P1/P2) ** (1.0d0/n_sw)
            T2_sw = T1 * (P2/P1) ** ((n_sw-1.0d0)/n_sw)
            W_sw = (P2*V2_sw - P1*V1) / (1.0d0 - n_sw)
        end if
        Q_sw = mass * cv * (T2_sw - T1) / 1000.0d0 + W_sw
        write(*,'(F8.4,4X,F12.2,4X,F12.4,4X,F12.4)') n_sw, T2_sw, W_sw, Q_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Polytropic relation: P*V^n = constant.'
    write(*,'(A)') '  V2 = V1*(P1/P2)^(1/n).'
    write(*,'(A)') '  T2 = T1*(P2/P1)^((n-1)/n).'
    write(*,'(A)') '  W = (P2V2 - P1V1)/(1-n) for n != 1.'
    write(*,'(A)') '  W = P1*V1*ln(V2/V1) for n = 1 (isothermal).'
    write(*,'(A)') '  Q = dU + W (first law, closed system).'
    write(*,'(A)') '  dS = m*[cv*ln(T2/T1) + R*ln(V2/V1)] (ideal gas).'
    write(*,'(A)') '  Special cases: n=0 isobaric, n=1 isothermal,'
    write(*,'(A)') '    n=gamma isentropic, n->inf isochoric.'

end program isentropic_polytropic
