program transonic_similarity
    implicit none
    double precision, parameter :: PI = 3.141592653589793d0
    double precision :: M_inf, gamma, tau, Cp_inc, CL_inc
    double precision :: chord, p_inf, T_inf, R_gas
    double precision :: gm1, gp1
    double precision :: rho_inf, a_inf, u_inf, q_inf
    double precision :: K_param, tau_23, beta_val, inv_beta
    double precision :: Cp_PG, Cp_bar, Cp_star, CL_PG
    double precision :: M_dd_conv, M_dd_super
    character(len=30) :: regime
    logical :: is_sonic
    integer :: i, n_points, iostat_val
    double precision :: dM, M_cur, K_cur, b_cur, ib_cur, cp_cur

    read(*,*,iostat=iostat_val) M_inf
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid Mach.'; stop; end if
    read(*,*,iostat=iostat_val) gamma
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid gamma.'; stop; end if
    read(*,*,iostat=iostat_val) tau
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid tau.'; stop; end if
    read(*,*,iostat=iostat_val) Cp_inc
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid Cp_inc.'; stop; end if
    read(*,*,iostat=iostat_val) CL_inc
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid CL_inc.'; stop; end if
    read(*,*,iostat=iostat_val) chord
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid chord.'; stop; end if
    read(*,*,iostat=iostat_val) p_inf
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid pressure.'; stop; end if
    read(*,*,iostat=iostat_val) T_inf
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid temperature.'; stop; end if
    read(*,*,iostat=iostat_val) R_gas
    if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid R_gas.'; stop; end if

    if (M_inf<=0.0d0) then; write(*,*) 'ERROR: Mach must be > 0.'; stop; end if
    if (gamma<=1.0d0) then; write(*,*) 'ERROR: gamma must be > 1.'; stop; end if
    if (tau<=0.0d0) then; write(*,*) 'ERROR: tau must be > 0.'; stop; end if
    if (chord<=0.0d0) then; write(*,*) 'ERROR: chord must be > 0.'; stop; end if
    if (p_inf<=0.0d0) then; write(*,*) 'ERROR: Pressure must be > 0.'; stop; end if
    if (T_inf<=0.0d0) then; write(*,*) 'ERROR: Temperature must be > 0.'; stop; end if
    if (R_gas<=0.0d0) then; write(*,*) 'ERROR: R_gas must be > 0.'; stop; end if

    gm1 = gamma - 1.0d0; gp1 = gamma + 1.0d0
    rho_inf = p_inf/(R_gas*T_inf)
    a_inf = sqrt(gamma*R_gas*T_inf)
    u_inf = M_inf*a_inf
    q_inf = 0.5d0*rho_inf*u_inf**2

    tau_23 = tau**(2.0d0/3.0d0)
    K_param = (1.0d0 - M_inf**2)/tau_23

    beta_val = sqrt(abs(1.0d0 - M_inf**2))
    if (beta_val < 1.0d-12) beta_val = 1.0d-12
    inv_beta = 1.0d0/beta_val

    if (abs(1.0d0 - M_inf**2) > 1.0d-6) then
        Cp_PG = Cp_inc/sqrt(abs(1.0d0 - M_inf**2))
    else
        Cp_PG = -9999.0d0
    end if
    Cp_bar = Cp_PG/tau_23

    if (M_inf < 1.0d0 .and. beta_val > 1.0d-6) then
        CL_PG = CL_inc/beta_val
    else
        CL_PG = CL_inc
    end if

    if (M_inf > 1.0d-6) then
        Cp_star = (2.0d0/(gamma*M_inf**2))* &
            (((2.0d0/gp1)*(1.0d0+0.5d0*gm1*M_inf**2))**(gamma/gm1) - 1.0d0)
    else
        Cp_star = -999.0d0
    end if

    is_sonic = .false.
    if (M_inf > 0.01d0 .and. M_inf < 1.0d0) then
        if (Cp_PG < Cp_star) is_sonic = .true.
    end if

    if (K_param > 2.0d0) then
        regime = 'SUBSONIC-LIKE'
    else if (K_param < -2.0d0) then
        regime = 'SUPERSONIC-LIKE'
    else
        regime = 'TRANSONIC'
    end if

    M_dd_conv  = 0.87d0 - CL_inc/10.0d0 - tau
    M_dd_super = 0.95d0 - CL_inc/10.0d0 - tau

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   TRANSONIC SIMILARITY PARAMETER CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- FREESTREAM CONDITIONS -----------------------------------'
    write(*,'(A,F12.6)')    '  Freestream Mach (M_inf) = ', M_inf
    write(*,'(A,F12.6)')    '  Specific Heat Ratio (g) = ', gamma
    write(*,'(A,ES14.6,A)') '  Static Pressure (p_inf) = ', p_inf, ' Pa'
    write(*,'(A,F12.2,A)')  '  Static Temperature      = ', T_inf, ' K'
    write(*,'(A,F12.4,A)')  '  Density (rho_inf)       = ', rho_inf, ' kg/m3'
    write(*,'(A,F12.4,A)')  '  Speed of Sound (a_inf)  = ', a_inf, ' m/s'
    write(*,'(A,F12.4,A)')  '  Velocity (u_inf)        = ', u_inf, ' m/s'
    write(*,'(A,ES14.6,A)') '  Dynamic Pressure (q)    = ', q_inf, ' Pa'
    write(*,*)
    write(*,'(A)') '--- AIRFOIL PARAMETERS --------------------------------------'
    write(*,'(A,F12.6)')    '  Thickness Ratio (tau)   = ', tau
    write(*,'(A,F12.8)')    '  tau^(2/3)               = ', tau_23
    write(*,'(A,F12.6)')    '  Cp (incompressible)     = ', Cp_inc
    write(*,'(A,F12.6)')    '  CL (incompressible)     = ', CL_inc
    write(*,'(A,F12.4,A)')  '  Chord Length            = ', chord, ' m'
    write(*,*)
    write(*,'(A)') '--- TRANSONIC SIMILARITY PARAMETER --------------------------'
    write(*,'(A,F14.8)')    '  K = (1-M^2)/tau^(2/3)   = ', K_param
    write(*,'(A,A)')        '  Flow Regime             = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- COMPRESSIBILITY CORRECTIONS -----------------------------'
    write(*,'(A,F12.8)')    '  beta = sqrt(|1-M^2|)    = ', beta_val
    write(*,'(A,F12.6)')    '  1/beta (amplification)  = ', inv_beta
    write(*,'(A,F14.8)')    '  Cp Prandtl-Glauert      = ', Cp_PG
    write(*,'(A,F14.8)')    '  Cp_bar = Cp_PG/tau^2/3  = ', Cp_bar
    write(*,'(A,F12.8)')    '  CL Prandtl-Glauert      = ', CL_PG
    write(*,*)
    write(*,'(A)') '--- CRITICAL CONDITION CHECK --------------------------------'
    write(*,'(A,F14.8)')    '  Cp* (sonic condition)   = ', Cp_star
    write(*,'(A,F14.8)')    '  Cp_PG at this Mach      = ', Cp_PG
    if (is_sonic) then
        write(*,'(A)')      '  Status: LOCALLY SUPERSONIC (Cp_PG < Cp*)'
    else
        write(*,'(A)')      '  Status: SUBSONIC EVERYWHERE'
    end if
    write(*,*)
    write(*,'(A)') '--- DRAG DIVERGENCE MACH (KORN EQUATION) --------------------'
    write(*,'(A,F12.6)')    '  M_dd (conventional)     = ', M_dd_conv
    write(*,'(A,F12.6)')    '  M_dd (supercritical)    = ', M_dd_super
    write(*,'(A)')          '  Korn: M_dd + CL/10 + t/c = kappa'
    write(*,'(A)')          '  kappa = 0.87 (conventional), 0.95 (supercritical)'
    if (M_inf > M_dd_conv) then
        write(*,'(A)')      '  *** DRAG DIVERGENCE: M_inf > M_dd(conv) ***'
    end if
    write(*,*)
    write(*,'(A)') '--- FORCE ESTIMATES (per unit span) -------------------------'
    write(*,'(A,ES14.6,A)') '  Lift / span             = ', q_inf*chord*CL_PG, ' N/m'
    write(*,*)

    write(*,'(A)') '--- PROFILE vs MACH NUMBER ----------------------------------'
    write(*,'(A)') '  M_inf       K             beta       1/beta     Cp_PG'
    write(*,'(A)') '  ----------------------------------------------------------'
    n_points = 40
    dM = 1.0d0/dble(n_points)
    do i = 0, n_points
        M_cur = 0.5d0 + dble(i)*dM
        K_cur = (1.0d0 - M_cur**2)/tau_23
        b_cur = sqrt(abs(1.0d0 - M_cur**2))
        if (b_cur < 1.0d-8) b_cur = 1.0d-8
        ib_cur = 1.0d0/b_cur
        if (abs(1.0d0 - M_cur**2) > 1.0d-6) then
            cp_cur = Cp_inc/sqrt(abs(1.0d0 - M_cur**2))
        else
            cp_cur = -9999.0d0
        end if
        write(*,'(F8.4,2X,F12.6,2X,F10.6,2X,F10.4,2X,F12.6)') &
            M_cur, K_cur, b_cur, ib_cur, cp_cur
    end do
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  K = (1 - M^2) / tau^(2/3)  (transonic similarity parameter)'
    write(*,'(A)') '  Cp_PG = Cp_0 / sqrt(|1 - M^2|)  (Prandtl-Glauert)'
    write(*,'(A)') '  Cp_bar = Cp / tau^(2/3)  (scaled pressure coefficient)'
    write(*,'(A)') '  Korn: M_dd + CL/10 + t/c = kappa'
    write(*,'(A)') '  TSD: (1-M^2)phi_xx + phi_yy = 0 (linear, breaks at M~1)'
    write(*,'(A)') '============================================================'
end program transonic_similarity