program ic_engine_cycle
    implicit none
    
    ! Constants for air
    double precision, parameter :: Cp = 1.005d0       ! kJ/(kg.K)
    double precision, parameter :: Cv = 0.718d0       ! kJ/(kg.K)
    double precision, parameter :: R = 0.287d0        ! kJ/(kg.K)
    double precision, parameter :: k_ratio = 1.4d0
    double precision, parameter :: T_ref = 273.15d0   ! K for entropy
    double precision, parameter :: P_ref = 100.0d0    ! kPa for entropy
    
    ! Inputs
    integer :: cycle_type  ! 1 = Otto, 2 = Diesel
    double precision :: T1_c, P1, r_ratio, T3_c, eta_c, eta_e, m_dot
    
    ! State properties (Temperatures in K, pressures in kPa, specific volumes in m3/kg, entropies in kJ/kgK)
    double precision :: T1, P1_val, v1, s1
    double precision :: T2, P2, v2, s2, T2s
    double precision :: T3, P3, v3, s3
    double precision :: T4, P4, v4, s4, T4s
    
    ! Cycle metrics
    double precision :: r_c        ! Cutoff ratio (Diesel only)
    double precision :: W_c, W_e, W_net
    double precision :: Q_in, Q_out, eta_thermal, MEP
    double precision :: W_net_total, Q_in_total
    
    ! Read from stdin
    read(*,*) cycle_type
    read(*,*) T1_c
    read(*,*) P1
    read(*,*) r_ratio
    read(*,*) T3_c
    read(*,*) eta_c
    read(*,*) eta_e
    read(*,*) m_dot
    
    ! Convert temperatures to Kelvin
    T1 = T1_c + 273.15d0
    T3 = T3_c + 273.15d0
    P1_val = P1
    
    ! -------------------------------------------------------------
    ! STATE 1: Compression Start
    ! -------------------------------------------------------------
    v1 = R * T1 / P1_val
    s1 = entropy(T1, P1_val)
    
    ! -------------------------------------------------------------
    ! STATE 2: Compression End
    ! -------------------------------------------------------------
    v2 = v1 / r_ratio
    T2s = T1 * (r_ratio)**(k_ratio - 1.0d0)
    T2 = T1 + (T2s - T1) / eta_c
    P2 = R * T2 / v2
    s2 = entropy(T2, P2)
    W_c = Cv * (T2 - T1)
    
    ! -------------------------------------------------------------
    ! STATE 3: Heat Addition End
    ! -------------------------------------------------------------
    if (cycle_type == 1) then
        ! Otto Cycle: Constant Volume Heat Addition
        v3 = v2
        P3 = P2 * (T3 / T2)
        s3 = entropy(T3, P3)
        Q_in = Cv * (T3 - T2)
        r_c = 1.0d0
    else
        ! Diesel Cycle: Constant Pressure Heat Addition
        P3 = P2
        r_c = T3 / T2  ! Cutoff ratio
        v3 = v2 * r_c
        s3 = entropy(T3, P3)
        Q_in = Cp * (T3 - T2)
    end if
    
    ! -------------------------------------------------------------
    ! STATE 4: Expansion End
    ! -------------------------------------------------------------
    v4 = v1
    if (cycle_type == 1) then
        ! Otto Cycle
        T4s = T3 / (r_ratio)**(k_ratio - 1.0d0)
        T4 = T3 - eta_e * (T3 - T4s)
        P4 = R * T4 / v4
        s4 = entropy(T4, P4)
        W_e = Cv * (T3 - T4)
        Q_out = Cv * (T4 - T1)
    else
        ! Diesel Cycle
        T4s = T3 / (r_ratio / r_c)**(k_ratio - 1.0d0)
        T4 = T3 - eta_e * (T3 - T4s)
        P4 = R * T4 / v4
        s4 = entropy(T4, P4)
        W_e = R * (T3 - T2) + Cv * (T3 - T4)
        Q_out = Cv * (T4 - T1)
    end if
    
    ! -------------------------------------------------------------
    ! PERFORMANCE METRICS
    ! -------------------------------------------------------------
    W_net = Q_in - Q_out
    eta_thermal = W_net / Q_in
    MEP = W_net / (v1 - v2)
    
    W_net_total = m_dot * W_net
    Q_in_total = m_dot * Q_in
    
    ! -------------------------------------------------------------
    ! OUTPUT IN PARSABLE FORMAT
    ! -------------------------------------------------------------
    write(*,'(A,I2)') 'CYCLE_TYPE=', cycle_type
    write(*,'(A,F14.4)') 'T1=', T1_c
    write(*,'(A,F14.4)') 'P1=', P1
    write(*,'(A,F14.6)') 'V1=', v1
    write(*,'(A,F14.6)') 'S1=', s1
    
    write(*,'(A,F14.4)') 'T2=', T2 - 273.15d0
    write(*,'(A,F14.4)') 'P2=', P2
    write(*,'(A,F14.6)') 'V2=', v2
    write(*,'(A,F14.6)') 'S2=', s2
    
    write(*,'(A,F14.4)') 'T3=', T3_c
    write(*,'(A,F14.4)') 'P3=', P3
    write(*,'(A,F14.6)') 'V3=', v3
    write(*,'(A,F14.6)') 'S3=', s3
    
    write(*,'(A,F14.4)') 'T4=', T4 - 273.15d0
    write(*,'(A,F14.4)') 'P4=', P4
    write(*,'(A,F14.6)') 'V4=', v4
    write(*,'(A,F14.6)') 'S4=', s4
    
    write(*,'(A,F14.4)') 'RC=', r_c
    write(*,'(A,F14.4)') 'WC=', W_c
    write(*,'(A,F14.4)') 'WE=', W_e
    write(*,'(A,F14.4)') 'WNET=', W_net
    write(*,'(A,F14.4)') 'QIN=', Q_in
    write(*,'(A,F14.4)') 'QOUT=', Q_out
    write(*,'(A,F14.4)') 'ETA_TH=', eta_thermal * 100.0d0
    write(*,'(A,F14.4)') 'MEP=', MEP
    
    write(*,'(A,F14.2)') 'WNET_TOT=', W_net_total
    write(*,'(A,F14.2)') 'QIN_TOT=', Q_in_total
    
contains

    ! Simple entropy calculation function
    double precision function entropy(T_k, P_kpa)
        double precision, intent(in) :: T_k, P_kpa
        entropy = Cp * log(T_k / T_ref) - R * log(P_kpa / P_ref)
    end function entropy

end program ic_engine_cycle