program brayton_cycle
    implicit none
    
    ! Constants for air (cold-air-standard assumptions)
    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 :: gamma = 1.4d0
    double precision, parameter :: T_ref = 273.15d0   ! K for entropy
    double precision, parameter :: P_ref = 100.0d0    ! kPa for entropy
    
    ! Inputs
    double precision :: T1_c, P1, r_p, T3_c, eta_c, eta_t, m_dot
    integer :: opt_regen, opt_inter, opt_reheat
    double precision :: eta_regen
    
    ! State properties (Temperatures in K, pressures in kPa, entropies in kJ/kgK)
    double precision :: T1, P1_val, s1
    double precision :: T1a, P1a, s1a, T1as
    double precision :: T1b, P1b, s1b
    double precision :: T2, P2, s2, T2s
    double precision :: T5, P5, s5
    double precision :: T3, P3, s3
    double precision :: T3a, P3a, s3a, T3as
    double precision :: T3b, P3b, s3b
    double precision :: T4, P4, s4, T4s
    double precision :: T6, P6, s6
    
    ! Cycle metrics
    double precision :: W_c1, W_c2, W_c, W_t1, W_t2, W_t, W_net
    double precision :: Q_in, Q_out, eta_thermal, BWR
    double precision :: W_net_total, Q_in_total
    
    ! Read from stdin
    read(*,*) T1_c
    read(*,*) P1
    read(*,*) r_p
    read(*,*) T3_c
    read(*,*) eta_c
    read(*,*) eta_t
    read(*,*) m_dot
    read(*,*) opt_regen
    read(*,*) eta_regen
    read(*,*) opt_inter
    read(*,*) opt_reheat
    
    ! Convert temperatures to Kelvin
    T1 = T1_c + 273.15d0
    T3 = T3_c + 273.15d0
    P1_val = P1
    
    ! -------------------------------------------------------------
    ! COMPRESSION PROCESS
    ! -------------------------------------------------------------
    s1 = entropy(T1, P1_val)
    
    if (opt_inter == 1) then
        ! 2-Stage Compression with Perfect Intercooling (intermediate ratio = sqrt(r_p))
        P1a = P1_val * sqrt(r_p)
        T1as = T1 * (P1a / P1_val)**((gamma - 1.0d0) / gamma)
        T1a = T1 + (T1as - T1) / eta_c
        s1a = entropy(T1a, P1a)
        W_c1 = Cp * (T1a - T1)
        
        ! Intercooling back to T1
        T1b = T1
        P1b = P1a
        s1b = entropy(T1b, P1b)
        
        ! Second compression stage
        P2 = P1_val * r_p
        T2s = T1b * (P2 / P1b)**((gamma - 1.0d0) / gamma)
        T2 = T1b + (T2s - T1b) / eta_c
        s2 = entropy(T2, P2)
        W_c2 = Cp * (T2 - T1b)
        
        W_c = W_c1 + W_c2
    else
        ! Single Stage Compression
        P2 = P1_val * r_p
        T2s = T1 * (P2 / P1_val)**((gamma - 1.0d0) / gamma)
        T2 = T1 + (T2s - T1) / eta_c
        s2 = entropy(T2, P2)
        W_c = Cp * (T2 - T1)
        
        ! Fill intermediate variables to avoid unitialized references
        T1a = T2
        P1a = P2
        s1a = s2
        T1b = T2
        P1b = P2
        s1b = s2
    end if
    
    ! -------------------------------------------------------------
    ! TURBINE EXPANSION PROCESS
    ! -------------------------------------------------------------
    P3 = P2
    s3 = entropy(T3, P3)
    
    if (opt_reheat == 1) then
        ! 2-Stage Expansion with Perfect Reheat
        P3a = P1_val * sqrt(r_p)
        T3as = T3 / (P3 / P3a)**((gamma - 1.0d0) / gamma)
        T3a = T3 - eta_t * (T3 - T3as)
        s3a = entropy(T3a, P3a)
        W_t1 = Cp * (T3 - T3a)
        
        ! Reheat back to T3
        T3b = T3
        P3b = P3a
        s3b = entropy(T3b, P3b)
        
        ! Second expansion stage
        P4 = P1_val
        T4s = T3b / (P3b / P4)**((gamma - 1.0d0) / gamma)
        T4 = T3b - eta_t * (T3b - T4s)
        s4 = entropy(T4, P4)
        W_t2 = Cp * (T3b - T4)
        
        W_t = W_t1 + W_t2
    else
        ! Single Stage Expansion
        P4 = P1_val
        T4s = T3 / (P3 / P4)**((gamma - 1.0d0) / gamma)
        T4 = T3 - eta_t * (T3 - T4s)
        s4 = entropy(T4, P4)
        W_t = Cp * (T3 - T4)
        
        ! Fill intermediate variables
        T3a = T4
        P3a = P4
        s3a = s4
        T3b = T4
        P3b = P4
        s3b = s4
    end if
    
    ! -------------------------------------------------------------
    ! REGENERATION
    ! -------------------------------------------------------------
    P5 = P2
    P6 = P4
    
    if (opt_regen == 1 .and. T4 > T2) then
        T5 = T2 + eta_regen * (T4 - T2)
        T6 = T4 - eta_regen * (T4 - T2)
    else
        T5 = T2
        T6 = T4
    end if
    
    s5 = entropy(T5, P5)
    s6 = entropy(T6, P6)
    
    ! -------------------------------------------------------------
    ! PERFORMANCE METRICS
    ! -------------------------------------------------------------
    W_net = W_t - W_c
    
    ! Heat input: Combustor (5 -> 3) + optional Reheater (3a -> 3b)
    Q_in = Cp * (T3 - T5)
    if (opt_reheat == 1) then
        Q_in = Q_in + Cp * (T3b - T3a)
    end if
    
    ! Heat rejected: Exhaust (6 -> 1) + optional Intercooler (1a -> 1b)
    Q_out = Cp * (T6 - T1)
    if (opt_inter == 1) then
        Q_out = Q_out + Cp * (T1a - T1b)
    end if
    
    eta_thermal = W_net / Q_in
    BWR = W_c / W_t
    
    W_net_total = m_dot * W_net
    Q_in_total = m_dot * Q_in
    
    ! -------------------------------------------------------------
    ! OUTPUT IN PARSABLE FORMAT
    ! -------------------------------------------------------------
    write(*,'(A,F14.4)') 'T1=', T1_c
    write(*,'(A,F14.4)') 'P1=', P1
    write(*,'(A,F14.4)') 'T2=', T2 - 273.15d0
    write(*,'(A,F14.4)') 'P2=', P2
    write(*,'(A,F14.4)') 'T3=', T3_c
    write(*,'(A,F14.4)') 'P3=', P3
    write(*,'(A,F14.4)') 'T4=', T4 - 273.15d0
    write(*,'(A,F14.4)') 'P4=', P4
    write(*,'(A,F14.4)') 'T5=', T5 - 273.15d0
    write(*,'(A,F14.4)') 'T6=', T6 - 273.15d0
    
    write(*,'(A,F14.4)') 'T1A=', T1a - 273.15d0
    write(*,'(A,F14.4)') 'P1A=', P1a
    write(*,'(A,F14.4)') 'T1B=', T1b - 273.15d0
    write(*,'(A,F14.4)') 'P1B=', P1b
    
    write(*,'(A,F14.4)') 'T3A=', T3a - 273.15d0
    write(*,'(A,F14.4)') 'P3A=', P3a
    write(*,'(A,F14.4)') 'T3B=', T3b - 273.15d0
    write(*,'(A,F14.4)') 'P3B=', P3b
    
    write(*,'(A,F14.6)') 'S1=', s1
    write(*,'(A,F14.6)') 'S2=', s2
    write(*,'(A,F14.6)') 'S3=', s3
    write(*,'(A,F14.6)') 'S4=', s4
    write(*,'(A,F14.6)') 'S5=', s5
    write(*,'(A,F14.6)') 'S6=', s6
    
    write(*,'(A,F14.4)') 'WC=', W_c
    write(*,'(A,F14.4)') 'WT=', W_t
    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)') 'BWR=', BWR * 100.0d0
    
    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 brayton_cycle
