program free_jet_impact
    implicit none
    integer :: i, iostat_val, n_points, vane_code
    double precision :: D0_mm, U0, rho, mu, X_mm, D0_m, X_m
    double precision :: A0, Q0, mdot0, Re0, F_thrust
    double precision :: L_core_m, L_core_mm, u_max_target, decay_factor
    double precision :: entrainment_ratio, Q_target, r_half_m, r_half_mm
    double precision :: q_dyn_Pa, P_stag_kPa, impact_force_N
    double precision :: x_pos_m, u_pos, r_curr_mm, eta, u_radial
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) D0_mm
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid nozzle diameter.'
        stop
    end if
    read(*,*,iostat=iostat_val) U0
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) X_mm
    read(*,*,iostat=iostat_val) vane_code ! 1 = flat 90 deg, 2 = inclined 45 deg, 3 = Pelton 180 deg

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all free jet inputs.'
        stop
    end if

    ! Validation
    if (D0_mm <= 0.0d0 .or. U0 <= 0.0d0 .or. rho <= 0.0d0 .or. mu <= 0.0d0 .or. X_mm <= 0.0d0) then
        write(*,*) 'ERROR: Geometry, velocity, and properties must be positive.'
        stop
    end if

    ! Conversions
    D0_m = D0_mm * 1.0d-3
    X_m = X_mm * 1.0d-3

    A0 = (PI / 4.0d0) * (D0_m**2)
    Q0 = A0 * U0
    mdot0 = rho * Q0
    Re0 = (rho * U0 * D0_m) / mu
    F_thrust = mdot0 * U0

    ! Potential Core
    L_core_m = 5.5d0 * D0_m
    L_core_mm = L_core_m * 1.0d3

    ! Centerline Velocity at Target Distance X
    if (X_m <= L_core_m) then
        u_max_target = U0
        entrainment_ratio = 1.0d0
    else
        decay_factor = 5.8d0 / (X_m / D0_m)
        u_max_target = min(U0, U0 * decay_factor)
        entrainment_ratio = max(1.0d0, 0.32d0 * (X_m / D0_m))
    end if

    Q_target = Q0 * entrainment_ratio
    r_half_m = 0.097d0 * X_m
    r_half_mm = max(D0_m / 2.0d0, r_half_m) * 1.0d3

    q_dyn_Pa = 0.5d0 * rho * (u_max_target**2)
    P_stag_kPa = q_dyn_Pa / 1.0d3

    ! Reaction Impact Force
    if (vane_code == 3) then
        ! Pelton 180 deg bucket
        impact_force_N = 2.0d0 * mdot0 * min(U0, u_max_target * 1.15d0) * 0.95d0
    else if (vane_code == 2) then
        ! 45 deg plate
        impact_force_N = mdot0 * min(U0, u_max_target * 1.15d0) * sin(45.0d0 * PI / 180.0d0)
    else
        ! Flat 90 deg plate
        impact_force_N = mdot0 * min(U0, u_max_target * 1.10d0)
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TURBULENT FREE JET & IMPINGEMENT ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,ES12.4,A)') 'Nozzle Exit Diameter D0   = ', D0_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Initial Jet Velocity U0   = ', U0, ' m/s'
    write(*,'(A,ES12.4)')   'Reynolds Number Re0       = ', Re0
    write(*,'(A,ES12.4,A)') 'Mass Flow Rate (mdot)     = ', mdot0, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Volumetric Flow Rate (Q0) = ', Q0, ' m3/s'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Potential Core Length     = ', L_core_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Target Distance X         = ', X_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Centerline Velocity u_max = ', u_max_target, ' m/s'
    write(*,'(A,ES12.4,A)') 'Jet Half-Width Radius r1/2= ', r_half_mm, ' mm'
    write(*,'(A,ES12.4)')   'Fluid Entrainment Q(X)/Q0 = ', entrainment_ratio
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Stagnation Pressure       = ', P_stag_kPa, ' kPa'
    write(*,'(A,ES12.4,A)') 'Nozzle Thrust Reaction F0 = ', F_thrust, ' N'
    write(*,'(A,ES12.4,A)') 'Target Impingement Force  = ', impact_force_N, ' N'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- AXIAL CENTERLINE VELOCITY DECAY U_MAX(X) ---'
    write(*,'(A)') '  x/D0         x_mm         u_max(m/s)'
    n_points = 20
    do i = 1, n_points
        x_pos_m = (D0_m * 0.5d0) + (35.0d0 * D0_m) * (dble(i - 1) / dble(n_points - 1))
        if (x_pos_m <= L_core_m) then
            u_pos = U0
        else
            u_pos = min(U0, U0 * (5.8d0 / (x_pos_m / D0_m)))
        end if
        write(*,'(2X,F10.2,2X,F10.2,2X,F10.2)') x_pos_m / D0_m, x_pos_m * 1.0d3, u_pos
    end do
    write(*,*)
    write(*,'(A)') '--- SCHLICHTING RADIAL VELOCITY PROFILE AT TARGET X ---'
    write(*,'(A)') '  r_mm         u(r)(m/s)'
    do i = -15, 15
        r_curr_mm = (r_half_mm * 2.2d0) * (dble(i) / 15.0d0)
        eta = abs(r_curr_mm) / max(0.1d0, r_half_mm)
        u_radial = u_max_target * ((1.0d0 + 0.414d0 * (eta**2))**(-2.0d0))
        write(*,'(2X,F10.2,2X,F10.2)') r_curr_mm, u_radial
    end do

end program free_jet_impact
