program prandtl_glauert
    implicit none
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: DEG2RAD = PI / 180.0d0
    double precision :: M_inf, gamma, Cp_inc, CL_inc, alpha_deg
    double precision :: chord, p_inf, T_inf, R_gas
    double precision :: gm1, gp1, alpha_rad
    double precision :: rho_inf, a_inf, u_inf, q_inf
    double precision :: beta, beta_sq, inv_beta
    double precision :: Cp_PG, Cp_KT, Cp_Laitone, denom_kt, denom_la
    double precision :: CL_PG, dCLda_PG
    double precision :: beta_sup, CL_Ack, Cd_wave, Cp_upper, Cp_lower
    double precision :: L_span, D_wave_span, Cp_star
    logical :: is_sub
    integer :: i, n_points, iostat_val
    double precision :: dM, M_cur, b_cur, cp_pg_val, cp_kt_val, cp_la_val, cl_cur, dk, dl

    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) 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) alpha_deg
    if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid alpha.'; 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 (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
    if (chord <= 0.0d0) then; write(*,*) 'ERROR: Chord must be > 0.'; stop; end if

    gm1 = gamma - 1.0d0
    gp1 = gamma + 1.0d0
    alpha_rad = alpha_deg * DEG2RAD
    is_sub = (M_inf < 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

    if (is_sub) then
        beta_sq = 1.0d0 - M_inf**2
        beta    = sqrt(beta_sq)
        inv_beta = 1.0d0 / beta
        Cp_PG = Cp_inc / beta
        CL_PG = CL_inc / beta
        dCLda_PG = 2.0d0 * PI / beta
        denom_kt = sqrt(beta_sq + M_inf**4 * ((Cp_inc + 1.0d0)**2) / 4.0d0)
        if (abs(denom_kt) > 1.0d-15) then
            Cp_KT = Cp_inc / denom_kt
        else
            Cp_KT = Cp_PG
        end if
        denom_la = beta + (M_inf**2 * (Cp_inc / 2.0d0)) / (1.0d0 + beta)
        if (abs(denom_la) > 1.0d-15) then
            Cp_Laitone = Cp_inc / denom_la
        else
            Cp_Laitone = Cp_PG
        end if
        L_span = q_inf * chord * CL_PG
        D_wave_span = 0.0d0
        Cp_star = (2.0d0 / (gamma * M_inf**2)) * &
            ( ((2.0d0/gp1)*(1.0d0+0.5d0*gm1*M_inf**2))**(gamma/gm1) - 1.0d0 )
        Cp_upper = 0.0d0; Cp_lower = 0.0d0; Cd_wave = 0.0d0; CL_Ack = 0.0d0; beta_sup = 0.0d0
    else
        beta_sup = sqrt(M_inf**2 - 1.0d0)
        beta = beta_sup; inv_beta = 1.0d0/beta
        Cp_upper = -2.0d0 * alpha_rad / beta_sup
        Cp_lower =  2.0d0 * alpha_rad / beta_sup
        Cp_PG = Cp_upper; Cp_KT = Cp_upper; Cp_Laitone = Cp_upper
        CL_Ack   = 4.0d0 * alpha_rad / beta_sup
        Cd_wave  = 4.0d0 * alpha_rad**2 / beta_sup
        CL_PG    = CL_Ack
        dCLda_PG = 4.0d0 / beta_sup
        L_span = q_inf * chord * CL_Ack
        D_wave_span = q_inf * chord * Cd_wave
        Cp_star = 0.0d0
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   PRANDTL-GLAUERT COMPRESSIBILITY CORRECTION 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)')  '  Freestream Velocity     = ', u_inf, ' m/s'
    write(*,'(A,ES14.6,A)') '  Dynamic Pressure (q)    = ', q_inf, ' Pa'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ----------------------------------------'
    write(*,'(A,F12.6)')    '  Cp (incompressible)     = ', Cp_inc
    write(*,'(A,F12.6)')    '  CL (incompressible)     = ', CL_inc
    write(*,'(A,F12.4,A)')  '  Angle of Attack (alpha) = ', alpha_deg, ' deg'
    write(*,'(A,F12.4,A)')  '  Chord Length            = ', chord, ' m'
    write(*,*)

    if (is_sub) then
        write(*,'(A)') '--- COMPRESSIBILITY FACTOR (SUBSONIC) -----------------------'
        write(*,'(A,F12.8)')    '  beta = sqrt(1-M^2)     = ', beta
        write(*,'(A,F12.8)')    '  1/beta (amplification)  = ', inv_beta
        write(*,*)
        write(*,'(A)') '--- PRESSURE COEFFICIENT CORRECTIONS ------------------------'
        write(*,'(A,F14.8)')    '  Cp Prandtl-Glauert     = ', Cp_PG
        write(*,'(A,F14.8)')    '  Cp Karman-Tsien         = ', Cp_KT
        write(*,'(A,F14.8)')    '  Cp Laitone              = ', Cp_Laitone
        write(*,'(A,F14.8)')    '  Cp* (sonic reference)   = ', Cp_star
        write(*,*)
        write(*,'(A)') '--- LIFT COEFFICIENT ----------------------------------------'
        write(*,'(A,F12.8)')    '  CL Prandtl-Glauert     = ', CL_PG
        write(*,'(A,F12.8)')    '  dCL/dalpha (per rad)    = ', dCLda_PG
        write(*,'(A,F12.8)')    '  dCL/dalpha (per deg)    = ', dCLda_PG * DEG2RAD
        write(*,*)
    else
        write(*,'(A)') '--- COMPRESSIBILITY FACTOR (SUPERSONIC) ---------------------'
        write(*,'(A,F12.8)')    '  beta = sqrt(M^2-1)     = ', beta_sup
        write(*,*)
        write(*,'(A)') '--- ACKERET LINEARIZED THEORY -------------------------------'
        write(*,'(A,F14.8)')    '  Cp upper surface        = ', Cp_upper
        write(*,'(A,F14.8)')    '  Cp lower surface        = ', Cp_lower
        write(*,'(A,F12.8)')    '  CL (Ackeret)            = ', CL_Ack
        write(*,'(A,F12.8)')    '  Cd_wave (Ackeret)       = ', Cd_wave
        write(*,'(A,F12.8)')    '  dCL/dalpha (per rad)    = ', dCLda_PG
        write(*,'(A,F12.8)')    '  dCL/dalpha (per deg)    = ', dCLda_PG * DEG2RAD
        write(*,*)
    end if

    write(*,'(A)') '--- FORCE ESTIMATES (per unit span) -------------------------'
    write(*,'(A,ES14.6,A)') '  Lift / span             = ', L_span, ' N/m'
    if (.not. is_sub) then
        write(*,'(A,ES14.6,A)') '  Wave Drag / span        = ', D_wave_span, ' N/m'
    end if
    write(*,*)

    write(*,'(A)') '--- CORRECTION PROFILE vs MACH NUMBER -----------------------'
    if (is_sub) then
        write(*,'(A)') '  M_inf     beta       Cp_PG       Cp_KT       Cp_Laitone  CL_PG'
    else
        write(*,'(A)') '  M_inf     beta_s     Cp_upper    Cp_lower    CL_Ack      Cd_wave'
    end if
    write(*,'(A)') '  ------------------------------------------------------------------'
    n_points = 40
    if (is_sub) then
        dM = 0.90d0 / dble(n_points)
        do i = 0, n_points
            M_cur = 0.05d0 + dble(i) * dM
            b_cur = sqrt(1.0d0 - M_cur**2)
            cp_pg_val = Cp_inc / b_cur
            dk = sqrt(1.0d0 - M_cur**2 + M_cur**4*((Cp_inc+1.0d0)**2)/4.0d0)
            if (abs(dk)>1.0d-15) then; cp_kt_val = Cp_inc/dk; else; cp_kt_val = cp_pg_val; end if
            dl = b_cur + (M_cur**2*(Cp_inc/2.0d0))/(1.0d0+b_cur)
            if (abs(dl)>1.0d-15) then; cp_la_val = Cp_inc/dl; else; cp_la_val = cp_pg_val; end if
            cl_cur = CL_inc / b_cur
            write(*,'(F8.4,2X,F10.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.6)') &
                M_cur, b_cur, cp_pg_val, cp_kt_val, cp_la_val, cl_cur
        end do
    else
        dM = (min(M_inf*1.5d0,5.0d0) - 1.05d0) / dble(n_points)
        do i = 0, n_points
            M_cur = 1.05d0 + dble(i) * dM
            b_cur = sqrt(M_cur**2 - 1.0d0)
            write(*,'(F8.4,2X,F10.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.8)') &
                M_cur, b_cur, &
                -2.0d0*alpha_rad/b_cur, 2.0d0*alpha_rad/b_cur, &
                4.0d0*alpha_rad/b_cur, 4.0d0*alpha_rad**2/b_cur
        end do
    end if
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Prandtl-Glauert: Cp(M) = Cp_0 / sqrt(1-M^2)'
    write(*,'(A)') '  Karman-Tsien:    Cp(M) = Cp_0 / sqrt(1-M^2+M^4*(Cp_0+1)^2/4)'
    write(*,'(A)') '  Laitone:         Cp(M) = Cp_0 / [beta + M^2*(Cp_0/2)/(1+beta)]'
    write(*,'(A)') '  Ackeret (super): Cp = +/-2*alpha/sqrt(M^2-1)'
    write(*,'(A)') '  CL_sub = CL_0/beta,  CL_sup = 4*alpha/sqrt(M^2-1)'
    write(*,'(A)') '============================================================'
end program prandtl_glauert
