program pitot_tube
    implicit none
    integer :: i, iostat_val
    double precision :: P_total, P_static, rho, gamma, R_gas, T_static
    double precision :: D_pipe, A_pipe, mu
    double precision :: V, q_dynamic, M, mdot, Re, P_impact
    double precision :: V_incomp, V_comp, correction, a_sound
    double precision :: dP_i, V_i, M_i, mdot_i, q_i
    double precision, parameter :: PI = 3.141592653589793d0
    character(len=60) :: comp_note

    read(*,*,iostat=iostat_val) P_total
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid total pressure input.'
        stop
    end if
    read(*,*,iostat=iostat_val) P_static
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) gamma
    read(*,*,iostat=iostat_val) R_gas
    read(*,*,iostat=iostat_val) T_static
    read(*,*,iostat=iostat_val) D_pipe
    read(*,*,iostat=iostat_val) mu
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all pitot tube inputs.'
        stop
    end if
    if (P_total <= 0.0d0 .or. P_static <= 0.0d0) then
        write(*,*) 'ERROR: Pressures must be positive.'
        stop
    end if
    if (P_total < P_static) then
        write(*,*) 'ERROR: Total pressure must be >= static pressure.'
        stop
    end if
    if (rho <= 0.0d0) then
        write(*,*) 'ERROR: Density must be positive.'
        stop
    end if
    if (gamma <= 1.0d0) gamma = 1.4d0
    if (R_gas <= 0.0d0) R_gas = 287.0d0
    if (T_static <= 0.0d0) T_static = 293.15d0
    if (D_pipe <= 0.0d0) D_pipe = 0.1d0
    if (mu <= 0.0d0) mu = 1.81d-5

    A_pipe = PI / 4.0d0 * D_pipe**2
    q_dynamic = P_total - P_static
    P_impact = q_dynamic

    ! Speed of sound
    a_sound = sqrt(gamma * R_gas * T_static)

    ! Incompressible Bernoulli: V = sqrt(2 dP / rho)
    V_incomp = sqrt(2.0d0 * q_dynamic / rho)

    ! Mach number from isentropic relation:
    ! P_total/P_static = (1 + (gamma-1)/2 M^2)^(gamma/(gamma-1))
    M = mach_from_pressure_ratio(P_total/P_static, gamma)

    ! Compressible velocity
    V_comp = M * a_sound

    ! Compressibility correction factor
    if (V_incomp > 1.0d-30) then
        correction = V_comp / V_incomp
    else
        correction = 1.0d0
    end if

    ! Choose appropriate velocity
    if (M < 0.3d0) then
        V = V_incomp
        comp_note = 'Incompressible (M < 0.3) — Bernoulli valid'
    else if (M < 1.0d0) then
        V = V_comp
        comp_note = 'Subsonic compressible — isentropic correction applied'
    else
        V = V_comp
        comp_note = 'Supersonic — Rayleigh pitot formula recommended'
    end if

    ! Mass flow rate
    mdot = rho * V * A_pipe

    ! Reynolds number
    Re = rho * V * D_pipe / mu

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   PITOT TUBE & PITOT-STATIC PROBE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Total Pressure Pt         = ', P_total, ' Pa'
    write(*,'(A,ES12.4,A)') '  Static Pressure Ps        = ', P_static, ' Pa'
    write(*,'(A,ES12.4,A)') '  Impact Pressure           = ', P_impact, ' Pa'
    write(*,'(A,ES12.4,A)') '  Fluid Density             = ', rho, ' kg/m3'
    write(*,'(A,ES12.4)')   '  Gamma                     = ', gamma
    write(*,'(A,ES12.4,A)') '  Gas Constant R            = ', R_gas, ' J/(kg.K)'
    write(*,'(A,ES12.4,A)') '  Static Temperature        = ', T_static, ' K'
    write(*,'(A,ES12.4,A)') '  Pipe Diameter             = ', D_pipe, ' m'
    write(*,'(A,ES12.4,A)') '  Viscosity                 = ', mu, ' Pa.s'
    write(*,*)
    write(*,'(A)') '--- VELOCITY RESULTS ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Speed of Sound            = ', a_sound, ' m/s'
    write(*,'(A,ES12.4,A)') '  Incompressible Velocity   = ', V_incomp, ' m/s'
    write(*,'(A,ES12.4,A)') '  Compressible Velocity     = ', V_comp, ' m/s'
    write(*,'(A,ES12.4)')   '  Compressibility Factor    = ', correction
    write(*,'(A,ES12.4,A)') '  Selected Velocity         = ', V, ' m/s'
    write(*,'(A,ES12.4)')   '  Mach Number               = ', M
    write(*,'(A,A)')        '  Flow Regime               = ', trim(comp_note)
    write(*,*)
    write(*,'(A)') '--- FLOW PARAMETERS -----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Dynamic Pressure q        = ', q_dynamic, ' Pa'
    write(*,'(A,ES12.4,A)') '  Mass Flow Rate            = ', mdot, ' kg/s'
    write(*,'(A,ES12.4)')   '  Reynolds Number           = ', Re
    write(*,'(A,ES12.4,A)') '  Pipe Area                 = ', A_pipe, ' m2'
    write(*,*)

    ! Velocity vs differential pressure sweep
    write(*,'(A)') '--- VELOCITY VS DELTA-P SWEEP -------------------------------'
    write(*,'(A)') '  dP[Pa]        V_incomp[m/s] V_comp[m/s]   Mach          correction'
    write(*,'(A)') '  -----------------------------------------------------------------------'
    do i = 1, 60
        dP_i = q_dynamic * 0.02d0 * dble(i)
        if (dP_i <= 0.0d0) cycle
        V_i = sqrt(2.0d0 * dP_i / rho)
        M_i = mach_from_pressure_ratio((P_static + dP_i)/P_static, gamma)
        if (M_i > 0.0d0) then
            write(*,'(ES12.4,2X,F12.4,2X,F12.4,2X,F10.5,2X,F10.5)') &
                dP_i, V_i, M_i*a_sound, M_i, &
                M_i*a_sound/max(V_i,1.0d-30)
        else
            write(*,'(ES12.4,2X,F12.4,2X,F12.4,2X,F10.5,2X,F10.5)') &
                dP_i, V_i, V_i, 0.0d0, 1.0d0
        end if
    end do
    write(*,*)

    ! Mach number vs pressure ratio sweep
    write(*,'(A)') '--- MACH VS PRESSURE RATIO SWEEP ----------------------------'
    write(*,'(A)') '  Pt/Ps         Mach          V[m/s]        T_total[K]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 50
        dP_i = 1.0d0 + 9.0d0 * dble(i-1) / 49.0d0   ! Pt/Ps from 1 to 10
        M_i = mach_from_pressure_ratio(dP_i, gamma)
        V_i = M_i * a_sound
        write(*,'(F10.4,2X,F10.5,2X,F12.3,2X,F12.2)') &
            dP_i, M_i, V_i, T_static*(1.0d0+(gamma-1.0d0)/2.0d0*M_i**2)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Bernoulli: V = sqrt(2 dP / rho).'
    write(*,'(A)') '  Isentropic: Pt/Ps = (1 + (gamma-1)/2 M^2)^(gamma/(gamma-1)).'
    write(*,'(A)') '  Compressible V = M * a; a = sqrt(gamma R T).'
    write(*,'(A)') '  Correction factor = V_comp / V_incomp.'

contains

    double precision function mach_from_pressure_ratio(pr, gam)
        implicit none
        double precision, intent(in) :: pr, gam
        double precision :: lo, hi, mid, pr_mid
        integer :: it
        if (pr <= 1.0d0) then
            mach_from_pressure_ratio = 0.0d0
            return
        end if
        lo = 0.0d0
        hi = 10.0d0
        do it = 1, 200
            mid = 0.5d0*(lo + hi)
            pr_mid = (1.0d0 + (gam-1.0d0)/2.0d0 * mid**2)**(gam/(gam-1.0d0))
            if (pr_mid < pr) then
                lo = mid
            else
                hi = mid
            end if
        end do
        mach_from_pressure_ratio = 0.5d0*(lo + hi)
    end function mach_from_pressure_ratio

end program pitot_tube