program compressible_flow
    implicit none
    
    double precision :: M, gamma, T0, P0, rho0
    double precision :: T_ratio, P_ratio, rho_ratio, A_ratio
    double precision :: a, V, T_s, P_s, rho_s, a_star
    double precision :: pi
    integer :: mode  ! 1=isentropic, 2=normal shock
    character(len=30) :: flow_type
    
    ! Normal shock variables
    double precision :: M2, P2_P1, T2_T1, rho2_rho1, P02_P01
    double precision :: dS_R
    
    integer :: i
    double precision :: M_i
    
    pi = 3.14159265358979d0
    
    ! Read inputs
    read(*,*) mode       ! 1=isentropic, 2=normal shock
    read(*,*) M          ! Mach number
    read(*,*) gamma      ! Specific heat ratio (1.4 for air)
    read(*,*) T0         ! Stagnation temperature [K]
    read(*,*) P0         ! Stagnation pressure [kPa]
    
    ! Stagnation density
    rho0 = P0 * 1000.0d0 / (287.0d0 * T0)  ! Assuming air (R=287 J/kg.K)
    
    ! ============================================
    ! ISENTROPIC FLOW RELATIONS
    ! ============================================
    T_ratio = 1.0d0 + (gamma - 1.0d0) / 2.0d0 * M**2
    T_s = T0 / T_ratio
    P_ratio = T_ratio**(gamma / (gamma - 1.0d0))
    P_s = P0 / P_ratio
    rho_ratio = T_ratio**(1.0d0 / (gamma - 1.0d0))
    rho_s = rho0 / rho_ratio
    
    ! Speed of sound
    a = sqrt(gamma * 287.0d0 * T_s)
    V = M * a
    a_star = sqrt(2.0d0 * gamma * 287.0d0 * T0 / (gamma + 1.0d0))
    
    ! Area ratio A/A*
    A_ratio = (1.0d0/M) * ((2.0d0/(gamma+1.0d0)) * (1.0d0 + (gamma-1.0d0)/2.0d0 * M**2)) &
              **((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
    
    if (M < 1.0d0) then
        flow_type = 'Subsonic'
    else if (abs(M - 1.0d0) < 1.0d-6) then
        flow_type = 'Sonic'
    else
        flow_type = 'Supersonic'
    end if
    
    ! ============================================
    ! OUTPUT — ISENTROPIC
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   COMPRESSIBLE FLOW CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- ISENTROPIC FLOW RELATIONS ---'
    write(*,'(A,F12.4)')    '  Mach Number (M)         = ', M
    write(*,'(A,F12.4)')    '  Specific Heat Ratio (gamma) = ', gamma
    write(*,'(A,A)')        '  Flow Type               = ', trim(flow_type)
    write(*,*)
    write(*,'(A)') '  Stagnation Conditions:'
    write(*,'(A,F12.2,A)')  '    T₀                    = ', T0, ' K'
    write(*,'(A,F12.4,A)')  '    P₀                    = ', P0, ' kPa'
    write(*,'(A,F12.4,A)')  '    rho₀                    = ', rho0, ' kg/m3'
    write(*,*)
    write(*,'(A)') '  Static Conditions:'
    write(*,'(A,F12.2,A)')  '    T                     = ', T_s, ' K'
    write(*,'(A,F12.4,A)')  '    P                     = ', P_s, ' kPa'
    write(*,'(A,F12.4,A)')  '    rho                     = ', rho_s, ' kg/m3'
    write(*,*)
    write(*,'(A)') '  Ratios:'
    write(*,'(A,F12.6)')    '    T/T₀                  = ', 1.0d0/T_ratio
    write(*,'(A,F12.6)')    '    P/P₀                  = ', 1.0d0/P_ratio
    write(*,'(A,F12.6)')    '    rho/rho₀                  = ', 1.0d0/rho_ratio
    write(*,'(A,F12.6)')    '    A/A*                  = ', A_ratio
    write(*,*)
    write(*,'(A)') '  Velocity & Speed of Sound:'
    write(*,'(A,F12.2,A)')  '    Speed of Sound (a)    = ', a, ' m/s'
    write(*,'(A,F12.2,A)')  '    Flow Velocity (V)     = ', V, ' m/s'
    write(*,'(A,F12.2,A)')  '    Critical a* (sonic)   = ', a_star, ' m/s'
    write(*,*)
    
    ! ============================================
    ! NORMAL SHOCK (if supersonic)
    ! ============================================
    if (mode == 2 .and. M > 1.0d0) then
        ! Normal shock relations
        M2 = sqrt(((gamma-1.0d0)*M**2 + 2.0d0) / (2.0d0*gamma*M**2 - (gamma-1.0d0)))
        P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M**2 - 1.0d0)
        T2_T1 = (1.0d0 + 2.0d0*gamma/(gamma+1.0d0)*(M**2-1.0d0)) * &
                ((2.0d0+(gamma-1.0d0)*M**2)/((gamma+1.0d0)*M**2))
        rho2_rho1 = ((gamma+1.0d0)*M**2) / ((gamma-1.0d0)*M**2 + 2.0d0)
        P02_P01 = (((gamma+1.0d0)*M**2)/((gamma-1.0d0)*M**2+2.0d0))**(gamma/(gamma-1.0d0)) &
                  * (1.0d0/(1.0d0+2.0d0*gamma/(gamma+1.0d0)*(M**2-1.0d0)))**(1.0d0/(gamma-1.0d0))
        
        ! Entropy change
        dS_R = -log(P02_P01)
        
        write(*,'(A)') '--- NORMAL SHOCK RELATIONS ---'
        write(*,'(A,F12.4)')    '  Upstream Mach (M1)     = ', M
        write(*,'(A,F12.6)')    '  Downstream Mach (M2)   = ', M2
        write(*,*)
        write(*,'(A)') '  Across Shock Ratios:'
        write(*,'(A,F12.6)')    '    P2/P1                = ', P2_P1
        write(*,'(A,F12.6)')    '    T2/T1                = ', T2_T1
        write(*,'(A,F12.6)')    '    rho2/rho1                = ', rho2_rho1
        write(*,'(A,F12.6)')    '    P₀2/P₀1              = ', P02_P01
        write(*,*)
        write(*,'(A)') '  Downstream Conditions:'
        write(*,'(A,F12.2,A)')  '    P2                   = ', P_s * P2_P1, ' kPa'
        write(*,'(A,F12.2,A)')  '    T2                   = ', T_s * T2_T1, ' K'
        write(*,'(A,F12.4)')    '    Entropy Change deltaS/R  = ', dS_R
        write(*,'(A,F12.2,A)')  '    Stag. Pressure Loss  = ', (1.0d0-P02_P01)*100, ' %'
        write(*,*)
    end if
    
    ! ============================================
    ! ISENTROPIC TABLE
    ! ============================================
    write(*,'(A)') '--- ISENTROPIC FLOW TABLE ---'
    write(*,'(A)') '  M         T/T₀        P/P₀        rho/rho₀        A/A*'
    write(*,'(A)') '  ---'
    
    do i = 1, 30
        M_i = 0.05d0 + dble(i-1) * 2.95d0 / 29.0d0
        
        T_ratio = 1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2
        
        write(*,'(F6.2,4X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6)') M_i, &
            1.0d0/T_ratio, &
            1.0d0/T_ratio**(gamma/(gamma-1.0d0)), &
            1.0d0/T_ratio**(1.0d0/(gamma-1.0d0)), &
            (1.0d0/M_i) * ((2.0d0/(gamma+1.0d0)) * T_ratio) &
              **((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
    end do
    
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  T₀/T = 1 + (gamma-1)/2 x m2'
    write(*,'(A)') '  P₀/P = (T₀/T)^(gamma/(gamma-1))'
    write(*,'(A)') '  rho₀/rho = (T₀/T)^(1/(gamma-1))'
    write(*,'(A)') '  A/A* = (1/M)[(2/(gamma+1))(1+(gamma-1)/2 m2)]^((gamma+1)/(2(gamma-1)))'
    
end program compressible_flow