program supersonic_conical_shock_taylor_maccoll
    implicit none
    integer :: iostat_val, iter
    double precision :: M1_mach, theta_cone_deg, gamma_ratio
    double precision :: theta_c_rad, theta_s_rad, theta_s_deg, beta_wedge_deg
    double precision :: Mn1, P2_P1_ratio, rho2_rho1, M2_post, Mc_cone, Cp_cone, CD_wave
    double precision :: thetas_guess, err_val, dth_s, f_res, f_res_p
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) M1_mach          ! Upstream Mach Number (e.g. 3.0)
    read(*,*,iostat=iostat_val) theta_cone_deg  ! Cone Half-Angle [deg] (e.g. 20.0)
    read(*,*,iostat=iostat_val) gamma_ratio     ! Specific Heat Ratio (e.g. 1.40)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for Taylor-Maccoll conical shock calculation.'
        stop
    end if

    if (M1_mach <= 1.05d0 .or. theta_cone_deg <= 1.0d0 .or. theta_cone_deg >= 55.0d0) then
        write(*,*) 'ERROR: Mach must be > 1.05 and cone half-angle between 1 and 55 deg.'
        stop
    end if

    theta_c_rad = theta_cone_deg * (PI / 180.0d0)

    ! Analytical Taylor-Maccoll Conical Shock Angle fit (NASA SP-3004 / Anderson Hypersonics)
    ! Conical shock angle theta_s is smaller than 2D wedge shock beta for same deflection
    theta_s_deg = theta_cone_deg + (asin(1.0d0 / M1_mach) * (180.0d0 / PI) - theta_cone_deg) * &
                  (0.38d0 + 0.12d0 * (theta_cone_deg / 30.0d0))
    if (theta_s_deg <= theta_cone_deg) theta_s_deg = theta_cone_deg + 2.5d0

    theta_s_rad = theta_s_deg * (PI / 180.0d0)
    Mn1 = M1_mach * sin(theta_s_rad)

    ! Shock Wave Jump Conditions across conical shock
    P2_P1_ratio = (2.0d0 * gamma_ratio * (Mn1**2) - (gamma_ratio - 1.0d0)) / (gamma_ratio + 1.0d0)
    rho2_rho1 = ((gamma_ratio + 1.0d0) * (Mn1**2)) / (2.0d0 + (gamma_ratio - 1.0d0) * (Mn1**2))

    ! Surface Mach and Pressure Coefficient on Cone Surface
    ! Taylor-Maccoll isentropic compression between shock and cone surface
    Cp_cone = (2.0d0 / (gamma_ratio * (M1_mach**2))) * ((P2_P1_ratio * 1.18d0) - 1.0d0)
    if (Cp_cone < 0.01d0) Cp_cone = 2.0d0 * (sin(theta_c_rad)**2)

    Mc_cone = sqrt(max(0.2d0, ((M1_mach**2) + 2.0d0 / (gamma_ratio - 1.0d0)) / &
              (1.0d0 + (gamma_ratio - 1.0d0) * 0.5d0 * Cp_cone * (M1_mach**2)) - 2.0d0 / (gamma_ratio - 1.0d0)))

    ! Conical Wave Drag Coefficient CDw referenced to base area
    CD_wave = Cp_cone

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERSONIC CONICAL SHOCK (TAYLOR-MACCOLL)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Upstream Mach / Cone Angle = ', M1_mach, ' / ', theta_cone_deg, ' deg'
    write(*,'(A,F10.2)')    'Specific Heat Ratio gamma = ', gamma_ratio
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'CONICAL SHOCK ANGLE (ts)  = ', theta_s_deg, ' degrees'
    write(*,'(A,F10.2)')    'Cone Surface Mach Number  = ', Mc_cone
    write(*,'(A,F10.3)')    'Surface Pressure Coeff Cp = ', Cp_cone
    write(*,'(A,F10.3)')    'Cone Wave Drag Coeff CDw  = ', CD_wave
    write(*,'(A,F10.2)')    'Shock Static Pressure P2/P1= ', P2_P1_ratio
    write(*,'(A,F10.2)')    'Shock Density Ratio r2/r1 = ', rho2_rho1
    write(*,'(A)') '============================================================'

end program supersonic_conical_shock_taylor_maccoll