Oblique Shock Wave relations

! ============================================================================
! ThermoFluidCalc — Oblique Shock Wave Relations Solver
! Reference: NACA Report 1135, Anderson Modern Compressible Flow Ch. 4
! ============================================================================
program oblique_shock
    implicit none
    
    ! Input variables
    integer :: solve_for ! 1=Solve for beta (Weak), 2=Solve for beta (Strong), 3=Solve for theta, 4=Solve for M1
    double precision :: M1 ! Upstream Mach number
    double precision :: theta ! Deflection angle [deg]
    double precision :: beta ! Shock angle [deg]
    double precision :: gamma ! Specific heat ratio (default 1.4)
    
    ! Constants
    double precision, parameter :: pi = 3.141592653589793d0
    double precision, parameter :: rad = pi / 180.0d0
    
    ! Intermediate and output variables
    double precision :: theta_rad, beta_rad, mu, mu_deg
    double precision :: beta_star_rad, beta_star_deg, theta_max_rad, theta_max_deg
    double precision :: Mn1, Mn2, M2
    double precision :: P2_P1, T2_T1, rho2_rho1, P02_P01
    logical :: detached
    
    ! Solver variables
    double precision :: low, high, mid, g_mid, val
    integer :: iter
    
    ! Read input from stdin
    read(*,*) solve_for
    read(*,*) M1
    read(*,*) theta
    read(*,*) beta
    read(*,*) gamma
    
    ! Check inputs
    if (gamma <= 1.0d0) gamma = 1.4d0
    detached = .false.
    
    ! Initialize defaults for output ratios
    P2_P1 = 1.0d0
    T2_T1 = 1.0d0
    rho2_rho1 = 1.0d0
    P02_P01 = 1.0d0
    M2 = 0.0d0
    Mn1 = 0.0d0
    Mn2 = 0.0d0
    mu = 0.0d0
    mu_deg = 0.0d0
    beta_star_rad = 0.0d0
    beta_star_deg = 0.0d0
    theta_max_deg = 0.0d0
    
    ! ── CASE 1, 2, 3: Upstream Mach checks ───────────────────
    if (solve_for /= 4) then
        if (M1 <= 1.0d0) then
            write(*,*) "ERROR: Upstream Mach number M1 must be supersonic (> 1.0) for oblique shocks."
            stop
        end if
        ! Calculate Mach Angle
        mu = asin(1.0d0 / M1)
        mu_deg = mu / rad
        
        ! Calculate beta* and theta_max analytically
        val = (gamma + 1.0d0) * (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2 + (gamma + 1.0d0)/16.0d0 * M1**4)
        val = (gamma + 1.0d0)/4.0d0 * M1**2 - 1.0d0 + sqrt(val)
        val = val / (gamma * M1**2)
        if (val >= 0.0d0 .and. val <= 1.0d0) then
            beta_star_rad = asin(sqrt(val))
            beta_star_deg = beta_star_rad / rad
            theta_max_rad = atan(f_beta(beta_star_rad, M1, gamma))
            theta_max_deg = theta_max_rad / rad
        end if
    end if
    
    ! ── SOLVE ACCORDING TO MODE ──────────────────────────────
    select case (solve_for)
        case (1, 2) ! Solve for beta (Weak or Strong)
            theta_rad = theta * rad
            if (theta > theta_max_deg) then
                detached = .true.
                beta_rad = beta_star_rad ! Default to max deflection limit
                beta = beta_star_deg
            else
                if (solve_for == 1) then
                    ! Weak solution on [mu, beta_star]
                    low = mu
                    high = beta_star_rad
                    do iter = 1, 150
                        mid = 0.5d0 * (low + high)
                        g_mid = atan(f_beta(mid, M1, gamma)) - theta_rad
                        if (g_mid < 0.0d0) then
                            low = mid
                        else
                            high = mid
                        end if
                        if (abs(high - low) < 1.0d-14) exit
                    end do
                    beta_rad = mid
                    beta = beta_rad / rad
                else
                    ! Strong solution on [beta_star, 90 deg]
                    low = beta_star_rad
                    high = 90.0d0 * rad
                    do iter = 1, 150
                        mid = 0.5d0 * (low + high)
                        g_mid = atan(f_beta(mid, M1, gamma)) - theta_rad
                        if (g_mid > 0.0d0) then
                            low = mid
                        else
                            high = mid
                        end if
                        if (abs(high - low) < 1.0d-14) exit
                    end do
                    beta_rad = mid
                    beta = beta_rad / rad
                end if
            end if
            
        case (3) ! Solve for theta given M1 and beta
            beta_rad = beta * rad
            if (beta_rad <= mu) then
                write(*,*) "ERROR: Shock angle beta must be greater than the Mach angle mu = ", mu_deg
                stop
            end if
            if (beta >= 90.0d0) then
                theta = 0.0d0
                theta_rad = 0.0d0
            else
                val = f_beta(beta_rad, M1, gamma)
                if (val < 0.0d0) then
                    write(*,*) "ERROR: Given shock angle results in non-physical negative deflection."
                    stop
                end if
                theta_rad = atan(val)
                theta = theta_rad / rad
            end if
            
        case (4) ! Solve for M1 given theta and beta
            theta_rad = theta * rad
            beta_rad = beta * rad
            if (beta <= theta) then
                write(*,*) "ERROR: Shock angle beta must be greater than deflection angle theta."
                stop
            end if
            
            ! Solve for M1^2 using rearranged equation
            val = sin(2.0d0 * beta_rad) - tan(theta_rad) * (gamma + cos(2.0d0 * beta_rad))
            if (val <= 0.0d0) then
                write(*,*) "ERROR: No physical Mach number solution exists for the given angles (denominator is non-positive)."
                stop
            end if
            
            mid = 2.0d0 * (1.0d0 / tan(beta_rad) + tan(theta_rad)) / val
            if (mid <= 1.0d0) then
                write(*,*) "ERROR: Calculated Mach number M1 must be supersonic (> 1.0)."
                stop
            end if
            
            M1 = sqrt(mid)
            mu = asin(1.0d0 / M1)
            mu_deg = mu / rad
            
            ! Calculate beta* and theta_max for solved M1
            val = (gamma + 1.0d0) * (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2 + (gamma + 1.0d0)/16.0d0 * M1**4)
            val = (gamma + 1.0d0)/4.0d0 * M1**2 - 1.0d0 + sqrt(val)
            val = val / (gamma * M1**2)
            if (val >= 0.0d0 .and. val <= 1.0d0) then
                beta_star_rad = asin(sqrt(val))
                beta_star_deg = beta_star_rad / rad
                theta_max_rad = atan(f_beta(beta_star_rad, M1, gamma))
                theta_max_deg = theta_max_rad / rad
            end if
            
        case default
            write(*,*) "ERROR: Invalid solve target mode."
            stop
    end select
    
    ! ── COMPUTE DOWNSTREAM PROPERTIES (IF ATTACHED) ──────────
    if (.not. detached) then
        beta_rad = beta * rad
        theta_rad = theta * rad
        Mn1 = M1 * sin(beta_rad)
        
        if (Mn1 > 1.0d0) then
            ! Normal Shock Relations applied to Mn1
            Mn2 = sqrt((2.0d0 + (gamma - 1.0d0) * Mn1**2) / (2.0d0 * gamma * Mn1**2 - (gamma - 1.0d0)))
            M2 = Mn2 / sin(beta_rad - theta_rad)
            
            P2_P1 = (2.0d0 * gamma * Mn1**2 - (gamma - 1.0d0)) / (gamma + 1.0d0)
            rho2_rho1 = ((gamma + 1.0d0) * Mn1**2) / ((gamma - 1.0d0) * Mn1**2 + 2.0d0)
            T2_T1 = P2_P1 / rho2_rho1
            
            P02_P01 = (rho2_rho1)**(gamma / (gamma - 1.0d0)) * (1.0d0 / P2_P1)**(1.0d0 / (gamma - 1.0d0))
        else
            ! Sonic/Subsonic normal component (should not occur if inputs are physical)
            M2 = M1
            P2_P1 = 1.0d0
            rho2_rho1 = 1.0d0
            T2_T1 = 1.0d0
            P02_P01 = 1.0d0
        end if
    end if
    
    ! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
    write(*, '(A, A)') "Detachment Status = ", merge("Detached", "Attached", detached)
    write(*, '(A, F14.6)') "Shock Angle (beta) = ", beta
    write(*, '(A, F14.6)') "Deflection Angle (theta) = ", theta
    write(*, '(A, F14.6)') "Upstream Mach (M1) = ", M1
    if (detached) then
        write(*, '(A, F14.6)') "Downstream Mach (M2) = ", 0.0d0
        write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", 0.0d0
        write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", 0.0d0
        write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", 0.0d0
        write(*, '(A, F14.6)') "Total Pressure Ratio (P02/P01) = ", 0.0d0
    else
        write(*, '(A, F14.6)') "Downstream Mach (M2) = ", M2
        write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
        write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
        write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
        write(*, '(A, F14.6)') "Total Pressure Ratio (P02/P01) = ", P02_P01
    end if
    write(*, '(A, F14.6)') "Max Deflection Angle (theta_max) = ", theta_max_deg
    write(*, '(A, F14.6)') "Shock Angle Max Deflection (beta_max) = ", beta_star_deg
    write(*, '(A, F14.6)') "Mach Angle (mu) = ", mu_deg
    write(*, '(A, F14.6)') "Upstream Normal Mach (Mn1) = ", Mn1
    write(*, '(A, F14.6)') "Downstream Normal Mach (Mn2) = ", Mn2
    write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma

contains

    ! Helper function to evaluate the theta-beta-Mach relation: tan(theta)
    double precision function f_beta(b, M, g)
        double precision, intent(in) :: b ! Shock angle [rad]
        double precision, intent(in) :: M ! Mach number
        double precision, intent(in) :: g ! Specific heat ratio
        
        f_beta = 2.0d0 * (1.0d0 / tan(b)) * (M**2 * (sin(b))**2 - 1.0d0) / &
                 (M**2 * (g + cos(2.0d0 * b)) + 2.0d0)
    end function f_beta

end program oblique_shock