program critical_cp
    implicit none

    double precision, parameter :: PI = 3.141592653589793d0
    double precision :: M_inf, gamma, Cp_min_incomp
    double precision :: gp1, gm1
    double precision :: Cp_star, Cp_vacuum, Cp_stag
    double precision :: p_star_ratio, T_star_ratio, rho_star_ratio, a_star_ratio
    double precision :: factor_inf, factor_sonic, p_ratio_sonic
    double precision :: M_cr, Cp_star_cr, Cp_min_cr
    logical :: has_Mcr
    double precision :: M_lo, M_hi, M_mid, f_lo, f_mid, f_hi
    integer :: iter, i, n_points
    integer, parameter :: MAX_ITER = 500
    double precision, parameter :: TOL = 1.0d-12
    double precision :: dM, M_cur, cps_cur, cpv_cur, cpst_cur, fac_cur, fac_s
    integer :: iostat_val

    read(*,*,iostat=iostat_val) M_inf
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Mach number.'; stop; end if
    read(*,*,iostat=iostat_val) gamma
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid gamma.'; stop; end if
    read(*,*,iostat=iostat_val) Cp_min_incomp
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Cp_min.'; 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

    gp1 = gamma + 1.0d0
    gm1 = gamma - 1.0d0

    factor_inf = 1.0d0 + 0.5d0 * gm1 * M_inf**2
    factor_sonic = 2.0d0 / gp1 * factor_inf
    p_ratio_sonic = factor_sonic ** (gamma / gm1)

    Cp_star = (2.0d0 / (gamma * M_inf**2)) * (p_ratio_sonic - 1.0d0)
    Cp_vacuum = -2.0d0 / (gamma * M_inf**2)
    Cp_stag = (2.0d0 / (gamma * M_inf**2)) * (factor_inf ** (gamma / gm1) - 1.0d0)

    T_star_ratio = factor_sonic
    p_star_ratio = p_ratio_sonic
    rho_star_ratio = p_star_ratio / T_star_ratio
    a_star_ratio = sqrt(T_star_ratio)

    has_Mcr = .false.
    M_cr = 0.0d0
    Cp_star_cr = 0.0d0
    Cp_min_cr = 0.0d0

    if (Cp_min_incomp < -1.0d-10) then
        M_lo = 0.01d0
        M_hi = 0.999d0
        f_lo = cp_star_func(M_lo, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_lo**2)
        f_hi = cp_star_func(M_hi, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_hi**2)
        if (f_lo * f_hi < 0.0d0) then
            do iter = 1, MAX_ITER
                M_mid = 0.5d0 * (M_lo + M_hi)
                f_mid = cp_star_func(M_mid, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_mid**2)
                if (abs(f_mid) < TOL .or. abs(M_hi - M_lo) < TOL) exit
                f_lo = cp_star_func(M_lo, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_lo**2)
                if (f_lo * f_mid < 0.0d0) then
                    M_hi = M_mid
                else
                    M_lo = M_mid
                end if
            end do
            M_cr = M_mid
            Cp_star_cr = cp_star_func(M_cr, gamma)
            Cp_min_cr = Cp_min_incomp / sqrt(1.0d0 - M_cr**2)
            has_Mcr = .true.
        end if
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CRITICAL PRESSURE COEFFICIENT CALCULATOR (Cp*)'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- FREESTREAM CONDITIONS -----------------------------------'
    write(*,'(A,F12.4)')       '  Freestream Mach (M_inf) = ', M_inf
    write(*,'(A,F12.6)')       '  Specific Heat Ratio (g) = ', gamma
    write(*,'(A,F12.6)')       '  Cp_min (incompressible) = ', Cp_min_incomp
    write(*,*)
    write(*,'(A)') '--- CRITICAL PRESSURE COEFFICIENT ---------------------------'
    write(*,'(A,F14.8)')       '  Cp*  (sonic condition)  = ', Cp_star
    write(*,'(A,F14.8)')       '  Cp_vacuum (p_local=0)   = ', Cp_vacuum
    write(*,'(A,F14.8)')       '  Cp_stagnation (M_loc=0) = ', Cp_stag
    write(*,*)
    write(*,'(A)') '--- SONIC CONDITION RATIOS ----------------------------------'
    write(*,'(A,F12.8)')       '  p*/p_inf                = ', p_star_ratio
    write(*,'(A,F12.8)')       '  T*/T_inf                = ', T_star_ratio
    write(*,'(A,F12.8)')       '  rho*/rho_inf            = ', rho_star_ratio
    write(*,'(A,F12.8)')       '  a*/a_inf                = ', a_star_ratio
    write(*,*)

    if (has_Mcr) then
        write(*,'(A)') '--- CRITICAL MACH NUMBER ------------------------------------'
        write(*,'(A,F12.8)')       '  Critical Mach M_cr      = ', M_cr
        write(*,'(A,F14.8)')       '  Cp* at M_cr             = ', Cp_star_cr
        write(*,'(A,F14.8)')       '  Cp_min at M_cr (P-G)    = ', Cp_min_cr
        write(*,*)
    else
        write(*,'(A)') '--- CRITICAL MACH NUMBER ------------------------------------'
        write(*,'(A)')             '  Critical Mach M_cr      = N/A (no intersection found)'
        write(*,*)
    end if

    write(*,'(A)') '--- Cp PROFILE vs FREESTREAM MACH ---------------------------'
    write(*,'(A)') '  M_inf       Cp*          Cp_vacuum    Cp_stag'
    write(*,'(A)') '  ----------------------------------------------------------'

    n_points = 40
    if (M_inf > 1.0d0) then
        dM = (min(M_inf * 1.5d0, 5.0d0) - 0.1d0) / dble(n_points)
    else
        dM = 0.95d0 / dble(n_points)
    end if

    do i = 0, n_points
        M_cur = 0.1d0 + dble(i) * dM
        if (M_cur < 0.01d0) M_cur = 0.01d0
        fac_cur = 1.0d0 + 0.5d0 * gm1 * M_cur**2
        fac_s = (2.0d0 / gp1 * fac_cur) ** (gamma / gm1)
        cps_cur = (2.0d0 / (gamma * M_cur**2)) * (fac_s - 1.0d0)
        cpv_cur = -2.0d0 / (gamma * M_cur**2)
        cpst_cur = (2.0d0 / (gamma * M_cur**2)) * (fac_cur ** (gamma / gm1) - 1.0d0)
        write(*,'(F8.4,2X,F12.8,2X,F12.8,2X,F12.8)') M_cur, cps_cur, cpv_cur, cpst_cur
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Cp* = (2/(g*M^2)) * [ (2/(g+1) * (1+(g-1)/2*M^2))^(g/(g-1)) - 1 ]'
    write(*,'(A)') '  Cp_vacuum = -2/(g*M^2)'
    write(*,'(A)') '  Cp_stag = (2/(g*M^2)) * [ (1+(g-1)/2*M^2)^(g/(g-1)) - 1 ]'
    write(*,'(A)') '  M_cr: solve Cp*(M) = Cp_min_inc / sqrt(1 - M^2)  (Prandtl-Glauert)'
    write(*,'(A)') '============================================================'

contains
    function cp_star_func(M, gam) result(cps)
        double precision, intent(in) :: M, gam
        double precision :: cps, fac
        fac = 1.0d0 + 0.5d0*(gam-1.0d0)*M**2
        cps = (2.0d0/(gam*M**2)) * ((2.0d0/(gam+1.0d0)*fac)**(gam/(gam-1.0d0)) - 1.0d0)
    end function cp_star_func

end program critical_cp
