!==============================================================================
! ThermoFluidCalc — Calculator #23 : Pressure Coefficient Cp
!==============================================================================
! Physics : Cp = (p - p_inf) / q_inf   where q_inf = 0.5 * rho * V_inf^2
!
! Modes:
!   1 = Single point evaluation
!   2 = Surface Cp distribution (N points)
!   3 = Compressible correction sweep (Prandtl-Glauert, Karman-Tsien, Laitone)
!
! Reference : Gupta, §5.7 and §7.8
!
! Build:
!   gfortran -O2 -o cp_pressure cp_pressure.f90
!==============================================================================
program cp_pressure
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: MAX_PTS = 5000

  integer  :: mode, N, npts, i
  real(dp) :: p, p_inf, rho, V_inf, gam, q_inf, Cp, Mach
  real(dp) :: xc(MAX_PTS), p_local(MAX_PTS), Cp_arr(MAX_PTS)
  real(dp) :: Cp_min, Cp_max, xc_cpmin, Cp_crit
  real(dp) :: Cp0, M_max, dM, M_cur
  real(dp) :: beta, Cp_PG, Cp_KT, Cp_Lait
  real(dp) :: term_kt, term_la
  character(len=40) :: interp

  read(*,*) mode

  select case (mode)

  !=========================================================================
  ! MODE 1 : Single point
  !=========================================================================
  case (1)
    backspace(5)
    read(*,*) mode, p, p_inf, rho, V_inf, gam

    if (rho <= 0.0_dp .or. V_inf <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Density and velocity must be positive.'
      stop
    end if

    q_inf = 0.5_dp * rho * V_inf**2
    Cp    = (p - p_inf) / q_inf
    Mach  = V_inf / sqrt(gam * p_inf / rho)

    ! Interpretation
    if (abs(Cp - 1.0_dp) < 0.05_dp) then
      interp = 'Stagnation point'
    else if (abs(Cp) < 0.01_dp) then
      interp = 'Freestream (neutral)'
    else if (Cp < 0.0_dp) then
      interp = 'Favorable pressure gradient (suction)'
    else
      interp = 'Adverse pressure gradient'
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Single Point'
    write(*,'(A,ES15.8)')   'P=', p
    write(*,'(A,ES15.8)')   'P_INF=', p_inf
    write(*,'(A,ES15.8)')   'RHO=', rho
    write(*,'(A,ES15.8)')   'V_INF=', V_inf
    write(*,'(A,ES15.8)')   'GAMMA=', gam
    write(*,'(A,ES15.8)')   'Q_INF=', q_inf
    write(*,'(A,F12.6)')    'CP=', Cp
    write(*,'(A,F12.6)')    'MACH=', Mach
    write(*,'(A,A)')        'INTERPRETATION=', trim(interp)

  !=========================================================================
  ! MODE 2 : Surface distribution
  !=========================================================================
  case (2)
    backspace(5)
    read(*,*) mode, p_inf, rho, V_inf, gam, N

    if (N < 1 .or. N > MAX_PTS) then
      write(*,'(A)') 'ERROR=Number of points must be 1-5000.'
      stop
    end if
    if (rho <= 0.0_dp .or. V_inf <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Density and velocity must be positive.'
      stop
    end if

    q_inf = 0.5_dp * rho * V_inf**2
    Mach  = V_inf / sqrt(gam * p_inf / rho)

    do i = 1, N
      read(*,*) xc(i), p_local(i)
      Cp_arr(i) = (p_local(i) - p_inf) / q_inf
    end do

    ! Find min and max Cp
    Cp_min    = Cp_arr(1)
    Cp_max    = Cp_arr(1)
    xc_cpmin  = xc(1)
    do i = 2, N
      if (Cp_arr(i) < Cp_min) then
        Cp_min   = Cp_arr(i)
        xc_cpmin = xc(i)
      end if
      if (Cp_arr(i) > Cp_max) Cp_max = Cp_arr(i)
    end do

    ! Critical Cp (isentropic)
    if (Mach > 0.0_dp .and. Mach < 1.0_dp) then
      Cp_crit = (2.0_dp / (gam * Mach**2)) * &
        ( ((2.0_dp + (gam - 1.0_dp)*Mach**2) / (gam + 1.0_dp))**(gam/(gam-1.0_dp)) &
          - 1.0_dp )
    else
      Cp_crit = -999.0_dp
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Surface Distribution'
    write(*,'(A,ES15.8)')   'Q_INF=', q_inf
    write(*,'(A,F12.6)')    'MACH=', Mach
    write(*,'(A,I5)')       'NPTS=', N
    write(*,'(A,F12.6)')    'CP_MIN=', Cp_min
    write(*,'(A,F12.8)')    'XC_AT_CPMIN=', xc_cpmin
    write(*,'(A,F12.6)')    'CP_MAX=', Cp_max
    write(*,'(A,F12.6)')    'CP_CRIT=', Cp_crit

    write(*,'(A)') 'DATA_START'
    do i = 1, N
      write(*,'(F12.8,A,F12.6)') xc(i), ',', Cp_arr(i)
    end do
    write(*,'(A)') 'DATA_END'

  !=========================================================================
  ! MODE 3 : Compressible correction sweep
  !=========================================================================
  case (3)
    backspace(5)
    read(*,*) mode, Cp0, gam, M_max, npts

    if (npts < 2) npts = 2
    if (npts > MAX_PTS) npts = MAX_PTS
    if (M_max >= 1.0_dp) M_max = 0.99_dp
    if (M_max <= 0.0_dp) M_max = 0.9_dp

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Compressible Correction'
    write(*,'(A,F12.6)')    'CP0=', Cp0
    write(*,'(A,F12.6)')    'GAMMA=', gam
    write(*,'(A,F12.6)')    'M_MAX=', M_max
    write(*,'(A,I5)')       'NPTS=', npts

    dM = M_max / real(npts - 1, dp)

    write(*,'(A)') 'DATA_START'
    do i = 0, npts - 1
      M_cur = real(i, dp) * dM
      if (M_cur >= 1.0_dp) M_cur = 0.999_dp

      beta = sqrt(max(1.0_dp - M_cur**2, 1.0e-12_dp))

      ! Prandtl-Glauert
      Cp_PG = Cp0 / beta

      ! Karman-Tsien
      term_kt = beta + (M_cur**2 * Cp0) / (2.0_dp * beta * (1.0_dp + beta))
      if (abs(term_kt) > 1.0e-15_dp) then
        Cp_KT = Cp0 / term_kt
      else
        Cp_KT = Cp0 / beta
      end if

      ! Laitone
      term_la = beta + (M_cur**2 * (1.0_dp + (gam - 1.0_dp)/2.0_dp * M_cur**2) * Cp0) &
                / (2.0_dp * beta)
      if (abs(term_la) > 1.0e-15_dp) then
        Cp_Lait = Cp0 / term_la
      else
        Cp_Lait = Cp0 / beta
      end if

      write(*,'(F10.6,3(A,F12.6))') M_cur, ',', Cp_PG, ',', Cp_KT, ',', Cp_Lait
    end do
    write(*,'(A)') 'DATA_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (must be 1-3).'
    stop
  end select

end program cp_pressure