! ============================================================================
! ThermoFluidCalc — Prandtl-Meyer Expansion Fan Solver
! Reference: Anderson Modern Compressible Flow Ch. 4
! ============================================================================
program prandtl_meyer
    implicit none
    
    ! Input variables
    integer :: solve_for ! 1 = Given M1 & theta, solve M2; 2 = Given M1 & M2, solve theta; 3 = Given M2 & theta, solve M1
    double precision :: M1 ! Upstream Mach number
    double precision :: theta ! Deflection angle [deg]
    double precision :: M2 ! Downstream Mach number
    double precision :: gamma ! Specific heat ratio (default 1.4)
    double precision :: T1 ! Upstream Temperature [K] (0 = skip)
    double precision :: P1 ! Upstream Pressure [kPa] (0 = skip)
    
    ! Constants
    double precision, parameter :: pi = 3.141592653589793d0
    double precision, parameter :: rad = pi / 180.0d0
    double precision, parameter :: R_air = 287.0d0 ! Gas constant for air [J/kg-K]
    
    ! Intermediate and output variables
    double precision :: nu1_rad, nu2_rad, theta_rad
    double precision :: mu1, mu2, nu_max_rad, nu_max_deg
    double precision :: P2_P1, T2_T1, rho2_rho1
    double precision :: T2, P2, a1, a2, V1, V2
    
    ! Read input from stdin
    read(*,*) solve_for
    read(*,*) M1
    read(*,*) theta
    read(*,*) M2
    read(*,*) gamma
    read(*,*) T1
    read(*,*) P1
    
    ! Check inputs
    if (gamma <= 1.0d0) gamma = 1.4d0
    
    ! Compute nu_max analytically for the given gamma
    nu_max_rad = (sqrt((gamma + 1.0d0) / (gamma - 1.0d0)) - 1.0d0) * (pi / 2.0d0)
    nu_max_deg = nu_max_rad / rad
    
    ! ── PERFORM SOLVING ACCORDING TO MODE ────────────────────
    select case (solve_for)
        case (1) ! Given M1 and theta, solve M2
            if (M1 <= 1.0d0) then
                write(*,*) "ERROR: Upstream Mach number M1 must be supersonic (> 1.0) for Prandtl-Meyer expansion."
                stop
            end if
            nu1_rad = nu_mach(M1, gamma)
            theta_rad = theta * rad
            nu2_rad = nu1_rad + theta_rad
            
            if (nu2_rad >= nu_max_rad) then
                write(*,*) "ERROR: Deflection angle theta exceeds maximum possible turn angle for M1 = ", M1
                stop
            end if
            
            M2 = solve_mach_from_nu(nu2_rad, gamma, M1)
            
        case (2) ! Given M1 and M2, solve theta
            if (M1 <= 1.0d0 .or. M2 <= 1.0d0) then
                write(*,*) "ERROR: Mach numbers must be supersonic (> 1.0) for expansion waves."
                stop
            end if
            if (M2 < M1) then
                write(*,*) "ERROR: Downstream Mach number M2 must be greater than upstream Mach number M1."
                stop
            end if
            
            nu1_rad = nu_mach(M1, gamma)
            nu2_rad = nu_mach(M2, gamma)
            theta_rad = nu2_rad - nu1_rad
            theta = theta_rad / rad
            
        case (3) ! Given M2 and theta, solve M1
            if (M2 <= 1.0d0) then
                write(*,*) "ERROR: Downstream Mach number M2 must be supersonic (> 1.0)."
                stop
            end if
            nu2_rad = nu_mach(M2, gamma)
            theta_rad = theta * rad
            nu1_rad = nu2_rad - theta_rad
            
            if (nu1_rad < 0.0d0) then
                write(*,*) "ERROR: Deflection angle theta exceeds total downstream PM angle."
                stop
            end if
            
            M1 = solve_mach_from_nu(nu1_rad, gamma, 1.5d0)
            
        case default
            write(*,*) "ERROR: Invalid solve target mode."
            stop
    end select
    
    ! ── COMPUTE DOWNSTREAM PROPERTIES ────────────────────────
    nu1_rad = nu_mach(M1, gamma)
    nu2_rad = nu_mach(M2, gamma)
    
    mu1 = asin(1.0d0 / M1) / rad
    mu2 = asin(1.0d0 / M2) / rad
    
    ! Isentropic ratios
    T2_T1 = (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2) / (1.0d0 + (gamma - 1.0d0)/2.0d0 * M2**2)
    P2_P1 = (T2_T1)**(gamma / (gamma - 1.0d0))
    rho2_rho1 = (T2_T1)**(1.0d0 / (gamma - 1.0d0))
    
    ! Optional static temperatures/pressures calculation
    T2 = 0.0d0
    P2 = 0.0d0
    V1 = 0.0d0
    V2 = 0.0d0
    if (T1 > 0.0d0) then
        T2 = T1 * T2_T1
        a1 = sqrt(gamma * R_air * T1)
        a2 = sqrt(gamma * R_air * T2)
        V1 = M1 * a1
        V2 = M2 * a2
    end if
    if (P1 > 0.0d0) then
        P2 = P1 * P2_P1
    end if
    
    ! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
    write(*, '(A, F14.6)') "Upstream Mach (M1) = ", M1
    write(*, '(A, F14.6)') "Downstream Mach (M2) = ", M2
    write(*, '(A, F14.6)') "Deflection Angle (theta) = ", theta
    write(*, '(A, F14.6)') "PM Angle Upstream (nu1) = ", nu1_rad / rad
    write(*, '(A, F14.6)') "PM Angle Downstream (nu2) = ", nu2_rad / rad
    write(*, '(A, F14.6)') "Mach Angle Upstream (mu1) = ", mu1
    write(*, '(A, F14.6)') "Mach Angle Downstream (mu2) = ", mu2
    write(*, '(A, F14.6)') "Max Turning Angle (nu_max) = ", nu_max_deg
    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)') "Specific Heat Ratio (gamma) = ", gamma
    
    if (T1 > 0.0d0) then
        write(*, '(A, F14.4)') "Upstream Temp (T1) = ", T1
        write(*, '(A, F14.4)') "Downstream Temp (T2) = ", T2
        write(*, '(A, F14.2)') "Upstream Velocity (V1) = ", V1
        write(*, '(A, F14.2)') "Downstream Velocity (V2) = ", V2
    end if
    if (P1 > 0.0d0) then
        write(*, '(A, F14.4)') "Upstream Pressure (P1) = ", P1
        write(*, '(A, F14.4)') "Downstream Pressure (P2) = ", P2
    end if

contains

    ! Prandtl-Meyer function calculation: nu(M) in radians
    double precision function nu_mach(M, g)
        double precision, intent(in) :: M, g
        double precision :: factor
        if (M <= 1.0d0) then
            nu_mach = 0.0d0
        else
            factor = sqrt((g + 1.0d0) / (g - 1.0d0))
            nu_mach = factor * atan(sqrt((g - 1.0d0) / (g + 1.0d0) * (M**2 - 1.0d0))) - &
                      atan(sqrt(M**2 - 1.0d0))
        end if
    end function nu_mach

    ! Derivative of the Prandtl-Meyer function: dnu/dM
    double precision function dnu_dM(M, g)
        double precision, intent(in) :: M, g
        if (M <= 1.0d0) then
            dnu_dM = 1.0d-10 ! Avoid division by zero
        else
            dnu_dM = sqrt(M**2 - 1.0d0) / (M * (1.0d0 + (g - 1.0d0) / 2.0d0 * M**2))
        end if
    end function dnu_dM

    ! Newton-Raphson solver to find M from a target nu in radians
    double precision function solve_mach_from_nu(nu_target, g, start_M)
        double precision, intent(in) :: nu_target, g, start_M
        double precision :: M, diff, f_val, df_val
        integer :: i
        
        if (nu_target <= 1.0d-12) then
            solve_mach_from_nu = 1.0d0
            return
        end if
        
        ! Initialize solver with starting Mach or standard guess
        M = start_M
        if (M <= 1.0d0) M = 2.0d0
        
        do i = 1, 100
            f_val = nu_mach(M, g) - nu_target
            df_val = dnu_dM(M, g)
            diff = f_val / df_val
            M = M - diff
            if (M < 1.0d0) M = 1.0001d0 ! Floor Mach to supersonic limit
            if (abs(diff) < 1.0d-14) exit
        end do
        
        solve_mach_from_nu = M
    end function solve_mach_from_nu

end program prandtl_meyer
