program mach_angle
    implicit none

    ! =================================================================
    !  Mach Angle / Mach Cone Calculator
    !  Reads: M, gamma, p1 [Pa], T1 [K], R_gas [J/kg.K], x_dist [m]
    ! =================================================================

    ! Constants
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: DEG2RAD = PI / 180.0d0
    double precision, parameter :: RAD2DEG = 180.0d0 / PI

    ! Inputs
    double precision :: M, gamma, p1, T1, R_gas, x_dist

    ! Derived quantities
    double precision :: rho1, a1, u1
    double precision :: gp1, gm1

    ! Mach angle
    double precision :: mu_rad, mu_deg, full_cone_deg
    double precision :: sin_mu, cos_mu, tan_mu

    ! Cone geometry at distance x
    double precision :: cone_radius, cone_diameter, zone_area

    ! Stagnation conditions
    double precision :: T0, p0, rho0

    ! Disturbance analysis
    double precision :: t_sound_1m, dist_source_1m, ratio_dist

    ! Profile sweep
    integer :: i, n_points
    double precision :: dM, M_cur, mu_cur, sin_cur, tan_cur, cos_cur
    double precision :: p_p0_cur, T_T0_cur, M_sweep_max

    integer :: iostat_val

    ! -----------------------------------------------------------------
    ! Read inputs
    ! -----------------------------------------------------------------
    read(*,*,iostat=iostat_val) M
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid Mach number input.'
        stop
    end if

    read(*,*,iostat=iostat_val) gamma
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid specific heat ratio input.'
        stop
    end if

    read(*,*,iostat=iostat_val) p1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid pressure input.'
        stop
    end if

    read(*,*,iostat=iostat_val) T1
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid temperature input.'
        stop
    end if

    read(*,*,iostat=iostat_val) R_gas
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid gas constant input.'
        stop
    end if

    read(*,*,iostat=iostat_val) x_dist
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid distance input.'
        stop
    end if

    ! -----------------------------------------------------------------
    ! Validate inputs
    ! -----------------------------------------------------------------
    if (M < 1.0d0) then
        write(*,*) 'ERROR: Mach number must be >= 1.0 for a Mach cone to exist.'
        stop
    end if
    if (gamma <= 1.0d0) then
        write(*,*) 'ERROR: Specific heat ratio gamma must be > 1.0.'
        stop
    end if
    if (p1 <= 0.0d0) then
        write(*,*) 'ERROR: Pressure must be positive.'
        stop
    end if
    if (T1 <= 0.0d0) then
        write(*,*) 'ERROR: Temperature must be positive.'
        stop
    end if
    if (R_gas <= 0.0d0) then
        write(*,*) 'ERROR: Gas constant must be positive.'
        stop
    end if
    if (x_dist < 0.0d0) then
        write(*,*) 'ERROR: Distance must be >= 0.'
        stop
    end if

    ! -----------------------------------------------------------------
    ! Convenience
    ! -----------------------------------------------------------------
    gp1 = gamma + 1.0d0
    gm1 = gamma - 1.0d0

    ! -----------------------------------------------------------------
    ! Flow state
    ! -----------------------------------------------------------------
    rho1 = p1 / (R_gas * T1)
    a1   = sqrt(gamma * R_gas * T1)
    u1   = M * a1

    ! -----------------------------------------------------------------
    ! Stagnation conditions
    ! -----------------------------------------------------------------
    T0   = T1 * (1.0d0 + 0.5d0 * gm1 * M**2)
    p0   = p1 * (1.0d0 + 0.5d0 * gm1 * M**2) ** (gamma / gm1)
    rho0 = rho1 * (1.0d0 + 0.5d0 * gm1 * M**2) ** (1.0d0 / gm1)

    ! -----------------------------------------------------------------
    ! Mach angle
    ! -----------------------------------------------------------------
    mu_rad = asin(1.0d0 / M)
    mu_deg = mu_rad * RAD2DEG
    full_cone_deg = 2.0d0 * mu_deg

    sin_mu = sin(mu_rad)
    cos_mu = cos(mu_rad)
    tan_mu = tan(mu_rad)

    ! -----------------------------------------------------------------
    ! Cone geometry at distance x
    ! -----------------------------------------------------------------
    if (x_dist > 0.0d0) then
        cone_radius   = x_dist * tan_mu
        cone_diameter = 2.0d0 * cone_radius
        zone_area     = PI * cone_radius**2
    else
        cone_radius   = 0.0d0
        cone_diameter = 0.0d0
        zone_area     = 0.0d0
    end if

    ! -----------------------------------------------------------------
    ! Disturbance analysis
    ! -----------------------------------------------------------------
    t_sound_1m   = 1.0d0 / a1                ! time for sound to travel 1 m
    dist_source_1m = u1 * t_sound_1m         ! distance source moves in that time
    ratio_dist   = a1 / u1                   ! = 1/M = sin(mu)

    ! -----------------------------------------------------------------
    ! Print results
    ! -----------------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   MACH ANGLE / MACH CONE CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- FLOW CONDITIONS -----------------------------------------'
    write(*,'(A,F12.4)')       '  Mach Number (M)         = ', M
    write(*,'(A,F12.6)')       '  Specific Heat Ratio (g) = ', gamma
    write(*,'(A,F12.2,A)')     '  Gas Constant (R)        = ', R_gas, ' J/kg.K'
    write(*,'(A,ES14.6,A)')    '  Static Pressure (p)     = ', p1, ' Pa'
    write(*,'(A,F12.2,A)')     '  Static Temperature (T)  = ', T1, ' K'
    write(*,'(A,F12.4,A)')     '  Density (rho)           = ', rho1, ' kg/m3'
    write(*,'(A,F12.4,A)')     '  Speed of Sound (a)      = ', a1, ' m/s'
    write(*,'(A,F12.4,A)')     '  Flow Velocity (u)       = ', u1, ' m/s'
    write(*,*)

    write(*,'(A)') '--- STAGNATION CONDITIONS -----------------------------------'
    write(*,'(A,ES14.6,A)')    '  Total Pressure p0       = ', p0, ' Pa'
    write(*,'(A,F12.2,A)')     '  Total Temperature T0    = ', T0, ' K'
    write(*,'(A,F12.4,A)')     '  Total Density rho0      = ', rho0, ' kg/m3'
    write(*,*)

    write(*,'(A)') '--- MACH ANGLE RESULTS --------------------------------------'
    write(*,'(A,F12.4,A)')     '  Mach Angle (mu)         = ', mu_deg, ' deg'
    write(*,'(A,F12.8,A)')     '  Mach Angle (mu)         = ', mu_rad, ' rad'
    write(*,'(A,F12.4,A)')     '  Full Cone Angle (2*mu)  = ', full_cone_deg, ' deg'
    write(*,'(A,F12.8)')       '  sin(mu) = 1/M           = ', sin_mu
    write(*,'(A,F12.8)')       '  cos(mu)                 = ', cos_mu
    write(*,'(A,F12.8)')       '  tan(mu)                 = ', tan_mu
    write(*,*)

    write(*,'(A)') '--- MACH CONE GEOMETRY --------------------------------------'
    write(*,'(A,F12.4,A)')     '  Reference Distance (x)  = ', x_dist, ' m'
    write(*,'(A,F12.6,A)')     '  Cone Radius at x        = ', cone_radius, ' m'
    write(*,'(A,F12.6,A)')     '  Cone Diameter at x      = ', cone_diameter, ' m'
    write(*,'(A,F12.6,A)')     '  Zone of Action Area     = ', zone_area, ' m2'
    write(*,*)

    write(*,'(A)') '--- DISTURBANCE PROPAGATION ---------------------------------'
    write(*,'(A,ES14.6,A)')    '  Time for sound (1 m)    = ', t_sound_1m, ' s'
    write(*,'(A,F12.4,A)')     '  Source travel (in t)    = ', dist_source_1m, ' m'
    write(*,'(A,F12.6)')       '  Ratio a/u = sin(mu)     = ', ratio_dist
    write(*,*)

    ! -----------------------------------------------------------------
    ! Profile sweep: M from 1.0 to M_sweep_max
    ! -----------------------------------------------------------------
    M_sweep_max = max(M * 2.0d0, 10.0d0)
    if (M_sweep_max < 2.0d0) M_sweep_max = 10.0d0

    write(*,'(A)') '--- MACH ANGLE PROFILE vs MACH NUMBER -----------------------'
    write(*,'(A)') '  M          mu [deg]    sin(mu)     tan(mu)' // &
                   '     cos(mu)     p/p0        T/T0'
    write(*,'(A)') '  ----------------------------------------------------------' // &
                   '----------------------------'

    n_points = 40
    dM = (M_sweep_max - 1.0d0) / dble(n_points)

    do i = 0, n_points
        M_cur = 1.0d0 + dble(i) * dM
        if (M_cur < 1.0d0) M_cur = 1.0d0

        mu_cur  = asin(1.0d0 / M_cur) * RAD2DEG
        sin_cur = 1.0d0 / M_cur
        cos_cur = sqrt(1.0d0 - sin_cur**2)
        if (cos_cur > 1.0d-15) then
            tan_cur = sin_cur / cos_cur
        else
            tan_cur = 0.0d0
        end if

        T_T0_cur = 1.0d0 / (1.0d0 + 0.5d0 * gm1 * M_cur**2)
        p_p0_cur = T_T0_cur ** (gamma / gm1)

        write(*,'(F8.4,2X,F11.4,2X,F11.8,2X,F11.8,2X,F11.8,2X,F11.8,2X,F11.8)') &
            M_cur, mu_cur, sin_cur, tan_cur, cos_cur, p_p0_cur, T_T0_cur
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  mu = asin(1/M)            (Mach angle)'
    write(*,'(A)') '  sin(mu) = 1/M = a/u       (wave-speed ratio)'
    write(*,'(A)') '  tan(mu) = 1/sqrt(M^2-1)   (cone radius/distance)'
    write(*,'(A)') '  Cone radius at x: r = x * tan(mu)'
    write(*,'(A)') '  Full cone angle: 2*mu'
    write(*,'(A)') '============================================================'

end program mach_angle
