program flow_cylinder
    implicit none
    integer :: i, iostat_val
    double precision :: D, Vinf, rho, mu, Re, St, f_vortex, CD, CL_osc
    double precision :: F_drag, F_lift_amp, omega_shed, V_lockin_lo, V_lockin_hi
    double precision :: fn_struct, zeta_struct, m_per_L
    double precision :: Re_i, CD_i, St_i, V_i, f_i
    double precision, parameter :: PI = 3.141592653589793d0
    character(len=60) :: regime

    read(*,*,iostat=iostat_val) D
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid cylinder diameter input.'
        stop
    end if
    read(*,*,iostat=iostat_val) Vinf
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) fn_struct
    read(*,*,iostat=iostat_val) zeta_struct
    read(*,*,iostat=iostat_val) m_per_L
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all cylinder flow inputs.'
        stop
    end if
    if (D <= 0.0d0 .or. Vinf <= 0.0d0) then
        write(*,*) 'ERROR: Diameter and velocity must be positive.'
        stop
    end if
    if (rho <= 0.0d0 .or. mu <= 0.0d0) then
        write(*,*) 'ERROR: Density and viscosity must be positive.'
        stop
    end if
    if (fn_struct <= 0.0d0) fn_struct = 1.0d0
    if (zeta_struct < 0.0d0) zeta_struct = 0.01d0
    if (m_per_L <= 0.0d0) m_per_L = 1.0d0

    Re = rho * Vinf * D / mu

    ! Strouhal number correlation
    St = strouhal_number(Re)

    ! Vortex shedding frequency
    f_vortex = St * Vinf / D
    omega_shed = 2.0d0 * PI * f_vortex

    ! Drag coefficient correlation
    CD = drag_coeff(Re)

    ! Oscillating lift coefficient amplitude (approximate)
    CL_osc = lift_osc_coeff(Re)

    ! Forces per unit length
    F_drag = 0.5d0 * rho * Vinf**2 * D * CD
    F_lift_amp = 0.5d0 * rho * Vinf**2 * D * CL_osc

    ! Lock-in velocity range (f_vortex ≈ fn_struct)
    ! Lock-in when 0.8 fn < f_shed < 1.2 fn  =>  V = f_n D / St
    if (St > 1.0d-8) then
        V_lockin_lo = 0.8d0 * fn_struct * D / St
        V_lockin_hi = 1.2d0 * fn_struct * D / St
    else
        V_lockin_lo = 0.0d0
        V_lockin_hi = 0.0d0
    end if

    ! Regime classification
    call classify_regime(Re, regime)

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   FLOW AROUND CYLINDERS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Cylinder Diameter         = ', D, ' m'
    write(*,'(A,ES12.4,A)') '  Freestream Velocity       = ', Vinf, ' m/s'
    write(*,'(A,ES12.4,A)') '  Fluid Density             = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity         = ', mu, ' Pa.s'
    write(*,'(A,ES12.4)')   '  Reynolds Number           = ', Re
    write(*,'(A,ES12.4,A)') '  Structural fn             = ', fn_struct, ' Hz'
    write(*,'(A,ES12.4)')   '  Structural Damping zeta   = ', zeta_struct
    write(*,'(A,ES12.4,A)') '  Mass per Length            = ', m_per_L, ' kg/m'
    write(*,*)
    write(*,'(A)') '--- VORTEX SHEDDING RESULTS ---------------------------------'
    write(*,'(A,ES12.4)')   '  Strouhal Number           = ', St
    write(*,'(A,ES12.4,A)') '  Shedding Frequency        = ', f_vortex, ' Hz'
    write(*,'(A,ES12.4,A)') '  Shedding Angular Freq     = ', omega_shed, ' rad/s'
    write(*,'(A,A)')        '  Reynolds Regime           = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- FORCE COEFFICIENTS --------------------------------------'
    write(*,'(A,ES12.4)')   '  Drag Coefficient CD       = ', CD
    write(*,'(A,ES12.4)')   '  Osc Lift Coefficient CL   = ', CL_osc
    write(*,'(A,ES12.4,A)') '  Drag Force per Length     = ', F_drag, ' N/m'
    write(*,'(A,ES12.4,A)') '  Lift Amplitude per Length = ', F_lift_amp, ' N/m'
    write(*,*)
    write(*,'(A)') '--- LOCK-IN ANALYSIS ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Lock-in Velocity Low      = ', V_lockin_lo, ' m/s'
    write(*,'(A,ES12.4,A)') '  Lock-in Velocity High     = ', V_lockin_hi, ' m/s'
    if (Vinf >= V_lockin_lo .and. Vinf <= V_lockin_hi) then
        write(*,'(A)') '  *** WARNING: CURRENT VELOCITY IS IN LOCK-IN RANGE ***'
    else
        write(*,'(A)') '  Current velocity is outside lock-in range.'
    end if
    write(*,*)

    ! CD, St vs Re sweep
    write(*,'(A)') '--- CD AND ST VS RE SWEEP -----------------------------------'
    write(*,'(A)') '  Re            CD            St            f_shed[Hz]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 80
        Re_i = 0.1d0 * (1.0d7/0.1d0)**(dble(i-1)/79.0d0)
        CD_i = drag_coeff(Re_i)
        St_i = strouhal_number(Re_i)
        V_i = Re_i * mu / (rho * D)
        f_i = St_i * V_i / D
        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') Re_i, CD_i, St_i, f_i
    end do
    write(*,*)

    ! Shedding frequency vs velocity sweep
    write(*,'(A)') '--- FREQUENCY VS VELOCITY SWEEP -----------------------------'
    write(*,'(A)') '  V[m/s]        f_shed[Hz]    f/fn          regime'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 50
        V_i = Vinf * 0.1d0 * dble(i)
        Re_i = rho * V_i * D / mu
        St_i = strouhal_number(Re_i)
        f_i = St_i * V_i / D
        if (f_i/fn_struct > 0.8d0 .and. f_i/fn_struct < 1.2d0) then
            write(*,'(ES12.4,2X,ES12.4,2X,F10.4,2X,A)') V_i, f_i, f_i/fn_struct, 'LOCK-IN'
        else
            write(*,'(ES12.4,2X,ES12.4,2X,F10.4,2X,A)') V_i, f_i, f_i/fn_struct, 'safe'
        end if
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  St = f D / V (Strouhal number).'
    write(*,'(A)') '  St ~ 0.198(1 - 19.7/Re) for subcritical Re.'
    write(*,'(A)') '  CD from piecewise empirical fit vs Re.'
    write(*,'(A)') '  Lock-in: 0.8 fn < f_shed < 1.2 fn.'

contains

    double precision function strouhal_number(Re_in)
        implicit none
        double precision, intent(in) :: Re_in
        if (Re_in < 40.0d0) then
            strouhal_number = 0.0d0   ! no periodic shedding
        else if (Re_in < 300.0d0) then
            strouhal_number = 0.198d0 * (1.0d0 - 19.7d0/Re_in)
        else if (Re_in < 3.0d5) then
            strouhal_number = 0.20d0   ! subcritical
        else if (Re_in < 3.5d6) then
            strouhal_number = 0.30d0   ! critical/supercritical
        else
            strouhal_number = 0.27d0   ! transcritical
        end if
        if (strouhal_number < 0.0d0) strouhal_number = 0.0d0
    end function strouhal_number

    double precision function drag_coeff(Re_in)
        implicit none
        double precision, intent(in) :: Re_in
        if (Re_in < 1.0d0) then
            ! Oseen-Stokes regime
            drag_coeff = 24.0d0 / max(Re_in, 1.0d-10) * &
                         (1.0d0 + 3.0d0/16.0d0*Re_in)
            if (drag_coeff > 200.0d0) drag_coeff = 200.0d0
        else if (Re_in < 1.0d3) then
            drag_coeff = 10.0d0 / sqrt(Re_in)
        else if (Re_in < 3.0d5) then
            drag_coeff = 1.2d0   ! subcritical
        else if (Re_in < 5.0d5) then
            ! drag crisis
            drag_coeff = 1.2d0 - 0.9d0*(Re_in-3.0d5)/2.0d5
        else if (Re_in < 3.5d6) then
            drag_coeff = 0.3d0   ! supercritical
        else
            drag_coeff = 0.6d0   ! transcritical recovery
        end if
        if (drag_coeff < 0.15d0) drag_coeff = 0.15d0
    end function drag_coeff

    double precision function lift_osc_coeff(Re_in)
        implicit none
        double precision, intent(in) :: Re_in
        if (Re_in < 40.0d0) then
            lift_osc_coeff = 0.0d0
        else if (Re_in < 300.0d0) then
            lift_osc_coeff = 0.3d0
        else if (Re_in < 3.0d5) then
            lift_osc_coeff = 0.5d0   ! subcritical
        else if (Re_in < 5.0d5) then
            lift_osc_coeff = 0.1d0   ! crisis suppression
        else
            lift_osc_coeff = 0.2d0   ! supercritical
        end if
    end function lift_osc_coeff

    subroutine classify_regime(Re_in, rname)
        implicit none
        double precision, intent(in) :: Re_in
        character(len=60), intent(out) :: rname
        if (Re_in < 5.0d0) then
            rname = 'Creeping flow — no separation'
        else if (Re_in < 40.0d0) then
            rname = 'Steady twin vortices — no shedding'
        else if (Re_in < 200.0d0) then
            rname = 'Laminar vortex street (2D periodic)'
        else if (Re_in < 300.0d0) then
            rname = 'Transition to 3D vortex shedding'
        else if (Re_in < 3.0d5) then
            rname = 'Subcritical — turbulent wake, laminar BL'
        else if (Re_in < 3.5d6) then
            rname = 'Critical / supercritical — drag crisis'
        else
            rname = 'Transcritical — turbulent BL and wake'
        end if
    end subroutine classify_regime

end program flow_cylinder