program rankine_hugoniot
    implicit none

    ! ---------------------------------------------------------------
    !  Normal Shock Relations  (Rankine-Hugoniot)
    !  Reads: M1, gamma, p1 [Pa], T1 [K], R_gas [J/kg.K]
    ! ---------------------------------------------------------------

    ! Inputs
    double precision :: M1, gamma, p1, T1, R_gas

    ! Derived upstream quantities
    double precision :: rho1, a1, u1

    ! Shock ratios
    double precision :: p_ratio, rho_ratio, T_ratio
    double precision :: M2_sq, M2
    double precision :: p0_ratio, ds_over_R, u_ratio

    ! Absolute downstream quantities
    double precision :: p2, T2, rho2, a2, u2

    ! Stagnation quantities
    double precision :: p01, p02, T01, T02

    ! Profile sweep variables
    integer :: i, n_points
    double precision :: dM, M_cur
    double precision :: pr_cur, rr_cur, Tr_cur, M2s_cur, M2_cur, p0r_cur, ds_cur

    ! Helpers
    double precision :: gp1, gm1, gp1h, gm1h   ! gamma+1, gamma-1, halves
    integer :: iostat_val

    ! ------------------------------------------------------------------
    ! Read inputs from stdin
    ! ------------------------------------------------------------------
    read(*,*,iostat=iostat_val) M1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid upstream Mach number input.'
        stop
    end if

    read(*,*,iostat=iostat_val) gamma
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid specific heat ratio input.'
        stop
    end if

    read(*,*,iostat=iostat_val) p1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid upstream pressure input.'
        stop
    end if

    read(*,*,iostat=iostat_val) T1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid upstream temperature input.'
        stop
    end if

    read(*,*,iostat=iostat_val) R_gas
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid gas constant input.'
        stop
    end if

    ! ------------------------------------------------------------------
    ! Validate inputs
    ! ------------------------------------------------------------------
    if (M1 < 1.0d0) then
        write(*,*) 'ERROR: Upstream Mach number must be >= 1.0 for a normal shock.'
        stop
    end if
    if (gamma <= 1.0d0) then
        write(*,*) 'ERROR: Specific heat ratio gamma must be > 1.0.'
        stop
    end if
    if (p1 <= 0.0d0) then
        write(*,*) 'ERROR: Upstream pressure must be positive.'
        stop
    end if
    if (T1 <= 0.0d0) then
        write(*,*) 'ERROR: Upstream temperature must be positive.'
        stop
    end if
    if (R_gas <= 0.0d0) then
        write(*,*) 'ERROR: Gas constant must be positive.'
        stop
    end if

    ! ------------------------------------------------------------------
    ! Convenience variables
    ! ------------------------------------------------------------------
    gp1  = gamma + 1.0d0
    gm1  = gamma - 1.0d0
    gp1h = gp1 / 2.0d0
    gm1h = gm1 / 2.0d0

    ! ------------------------------------------------------------------
    ! Upstream state
    ! ------------------------------------------------------------------
    rho1 = p1 / (R_gas * T1)
    a1   = sqrt(gamma * R_gas * T1)
    u1   = M1 * a1

    ! ------------------------------------------------------------------
    ! Normal shock ratios
    ! ------------------------------------------------------------------
    ! Static pressure ratio  p2/p1
    p_ratio = 1.0d0 + (2.0d0 * gamma / gp1) * (M1**2 - 1.0d0)

    ! Density ratio  rho2/rho1  (= u1/u2 from continuity)
    rho_ratio = (gp1 * M1**2) / (gm1 * M1**2 + 2.0d0)

    ! Temperature ratio  T2/T1
    T_ratio = p_ratio / rho_ratio

    ! Downstream Mach number squared
    M2_sq = (gm1 * M1**2 + 2.0d0) / (2.0d0 * gamma * M1**2 - gm1)
    M2    = sqrt(M2_sq)

    ! Velocity ratio  u2/u1  (inverse of density ratio from continuity)
    u_ratio = 1.0d0 / rho_ratio

    ! Total pressure ratio  p02/p01
    p0_ratio = (rho_ratio ** (gamma / gm1)) &
             * (p_ratio ** (-1.0d0 / gm1))

    ! Entropy change  Delta_s / R
    ds_over_R = -log(p0_ratio)

    ! ------------------------------------------------------------------
    ! Absolute downstream state
    ! ------------------------------------------------------------------
    p2   = p_ratio * p1
    T2   = T_ratio * T1
    rho2 = rho_ratio * rho1
    a2   = sqrt(gamma * R_gas * T2)
    u2   = M2 * a2

    ! ------------------------------------------------------------------
    ! Stagnation (total) conditions
    ! ------------------------------------------------------------------
    p01  = p1 * (1.0d0 + gm1h * M1**2) ** (gamma / gm1)
    T01  = T1 * (1.0d0 + gm1h * M1**2)
    p02  = p0_ratio * p01
    T02  = T01   ! Total temperature is constant across normal shock

    ! ------------------------------------------------------------------
    ! Print results
    ! ------------------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   NORMAL SHOCK RELATIONS  (Rankine-Hugoniot)'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- UPSTREAM CONDITIONS -------------------------------------'
    write(*,'(A,F12.4)')       '  Upstream Mach (M1)      = ', M1
    write(*,'(A,F12.6)')       '  Specific Heat Ratio (g) = ', gamma
    write(*,'(A,F12.2,A)')     '  Gas Constant (R)        = ', R_gas, ' J/kg.K'
    write(*,'(A,ES14.6,A)')    '  Upstream Pressure (p1)  = ', p1, ' Pa'
    write(*,'(A,F12.2,A)')     '  Upstream Temperature(T1)= ', T1, ' K'
    write(*,'(A,F12.4,A)')     '  Upstream Density (rho1) = ', rho1, ' kg/m3'
    write(*,'(A,F12.4,A)')     '  Speed of Sound (a1)     = ', a1, ' m/s'
    write(*,'(A,F12.4,A)')     '  Flow Velocity (u1)      = ', u1, ' m/s'
    write(*,*)

    write(*,'(A)') '--- STAGNATION (TOTAL) CONDITIONS ---------------------------'
    write(*,'(A,ES14.6,A)')    '  Total Pressure p01      = ', p01, ' Pa'
    write(*,'(A,F12.2,A)')     '  Total Temperature T01   = ', T01, ' K'
    write(*,*)

    write(*,'(A)') '--- NORMAL SHOCK RATIOS -------------------------------------'
    write(*,'(A,F12.6)')       '  p2/p1  (pressure)       = ', p_ratio
    write(*,'(A,F12.6)')       '  rho2/rho1 (density)     = ', rho_ratio
    write(*,'(A,F12.6)')       '  T2/T1  (temperature)    = ', T_ratio
    write(*,'(A,F12.6)')       '  u2/u1  (velocity)       = ', u_ratio
    write(*,'(A,F12.6)')       '  M2     (downstream Mach)= ', M2
    write(*,'(A,F12.8)')       '  p02/p01 (total pressure)= ', p0_ratio
    write(*,'(A,F12.8)')       '  Ds/R   (entropy rise)   = ', ds_over_R
    write(*,*)

    write(*,'(A)') '--- DOWNSTREAM CONDITIONS -----------------------------------'
    write(*,'(A,ES14.6,A)')    '  Downstream Pressure p2  = ', p2, ' Pa'
    write(*,'(A,F12.2,A)')     '  Downstream Temp T2      = ', T2, ' K'
    write(*,'(A,F12.4,A)')     '  Downstream Density rho2 = ', rho2, ' kg/m3'
    write(*,'(A,F12.4,A)')     '  Speed of Sound a2       = ', a2, ' m/s'
    write(*,'(A,F12.4,A)')     '  Flow Velocity u2        = ', u2, ' m/s'
    write(*,'(A,F12.6)')       '  Downstream Mach M2      = ', M2
    write(*,'(A,ES14.6,A)')    '  Total Pressure p02      = ', p02, ' Pa'
    write(*,'(A,F12.2,A)')     '  Total Temperature T02   = ', T02, ' K'
    write(*,*)

    ! ------------------------------------------------------------------
    ! Profile sweep:  M1 from 1.0 to user M1 in 40 steps
    ! ------------------------------------------------------------------
    write(*,'(A)') '--- SHOCK PROPERTY PROFILE vs MACH NUMBER -------------------'
    write(*,'(A)') '  M1         p2/p1       rho2/rho1   T2/T1       M2' // &
                   '          p02/p01     Ds/R'
    write(*,'(A)') '  ----------------------------------------------------------' // &
                   '----------------------------'

    n_points = 40
    if (M1 <= 1.0d0) then
        dM = 0.0d0
    else
        dM = (M1 - 1.0d0) / dble(n_points)
    end if

    do i = 0, n_points
        M_cur = 1.0d0 + dble(i) * dM
        if (M_cur < 1.0d0) M_cur = 1.0d0

        pr_cur  = 1.0d0 + (2.0d0 * gamma / gp1) * (M_cur**2 - 1.0d0)
        rr_cur  = (gp1 * M_cur**2) / (gm1 * M_cur**2 + 2.0d0)
        Tr_cur  = pr_cur / rr_cur
        M2s_cur = (gm1 * M_cur**2 + 2.0d0) / (2.0d0 * gamma * M_cur**2 - gm1)
        M2_cur  = sqrt(M2s_cur)
        p0r_cur = (rr_cur ** (gamma / gm1)) * (pr_cur ** (-1.0d0 / gm1))
        ds_cur  = -log(p0r_cur)

        write(*,'(F8.4,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.8,2X,F11.8)') &
            M_cur, pr_cur, rr_cur, Tr_cur, M2_cur, p0r_cur, ds_cur
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  p2/p1     = 1 + 2*gamma/(gamma+1) * (M1^2 - 1)'
    write(*,'(A)') '  rho2/rho1 = (gamma+1)*M1^2 / ((gamma-1)*M1^2 + 2)'
    write(*,'(A)') '  T2/T1     = (p2/p1) / (rho2/rho1)'
    write(*,'(A)') '  M2^2      = ((gamma-1)*M1^2 + 2) / (2*gamma*M1^2 - (gamma-1))'
    write(*,'(A)') '  p02/p01   = (rho2/rho1)^(g/(g-1)) * (p2/p1)^(-1/(g-1))'
    write(*,'(A)') '  Ds/R      = -ln(p02/p01)'
    write(*,'(A)') '============================================================'

end program rankine_hugoniot
