program drag_force_calculator
    implicit none
    
    double precision :: U_inf, rho, mu, nu, A_ref, D_char
    double precision :: Re, Cd, F_drag, q_dynamic
    double precision :: pi
    integer :: shape  ! 1=sphere, 2=cylinder, 3=flat-plate-normal, 4=streamlined
    character(len=30) :: shape_name
    integer :: i
    double precision :: Re_i, Cd_i
    
    pi = 3.14159265358979d0
    
    ! Read inputs
    read(*,*) shape      ! Shape type
    read(*,*) D_char     ! Characteristic dimension [m] (diameter or length)
    read(*,*) U_inf      ! Free stream velocity [m/s]
    read(*,*) rho        ! Fluid density [kg/m3]
    read(*,*) mu         ! Dynamic viscosity [Pa.s]
    
    nu = mu / rho
    Re = rho * U_inf * D_char / mu
    q_dynamic = 0.5d0 * rho * U_inf**2
    
    ! ============================================
    ! DRAG COEFFICIENT CORRELATIONS
    ! ============================================
    select case(shape)
    case(1)
        shape_name = 'Sphere'
        A_ref = pi * D_char**2 / 4.0d0
        
        if (Re < 1.0d0) then
            ! Stokes flow
            Cd = 24.0d0 / Re
        else if (Re < 1000.0d0) then
            ! Schiller-Naumann
            Cd = 24.0d0/Re * (1.0d0 + 0.15d0 * Re**0.687d0)
        else if (Re < 2.0d5) then
            ! Standard drag
            Cd = 0.44d0
        else
            ! Post-critical (turbulent boundary layer)
            Cd = 0.1d0
        end if
        
    case(2)
        shape_name = 'Infinite Cylinder'
        A_ref = D_char * 1.0d0  ! Per unit length
        
        if (Re < 1.0d0) then
            Cd = 8.0d0 * pi / (Re * (2.002d0 - log(Re)))
        else if (Re < 1000.0d0) then
            Cd = 1.0d0 + 10.0d0 * Re**(-2.0d0/3.0d0)
        else if (Re < 2.0d5) then
            Cd = 1.2d0
        else
            Cd = 0.3d0
        end if
        
    case(3)
        shape_name = 'Flat Plate (normal)'
        A_ref = D_char * D_char  ! Square plate
        Cd = 1.17d0  ! Approximately constant for Re > 1000
        if (Re < 100.0d0) then
            Cd = 2.0d0  ! Low Re approximation
        end if
        
    case(4)
        shape_name = 'Streamlined Body'
        A_ref = pi * (D_char/2.0d0)**2  ! Frontal area
        Cd = 0.04d0  ! Well-streamlined body
        if (Re < 1.0d5) then
            Cd = 0.1d0
        end if
        
    case default
        shape_name = 'Sphere'
        A_ref = pi * D_char**2 / 4.0d0
        Cd = 0.44d0
    end select
    
    ! Drag force
    F_drag = Cd * q_dynamic * A_ref
    
    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   EXTERNAL FLOW — DRAG FORCE CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- CONFIGURATION ---'
    write(*,'(A,A)')        '  Body Shape              = ', trim(shape_name)
    write(*,'(A,F12.4,A)')  '  Characteristic Dim (D)  = ', D_char, ' m'
    write(*,'(A,ES12.4,A)') '  Reference Area          = ', A_ref, ' m2'
    write(*,*)
    write(*,'(A)') '--- FLOW CONDITIONS ---'
    write(*,'(A,F12.4,A)')  '  Free Stream Velocity    = ', U_inf, ' m/s'
    write(*,'(A,F12.2,A)')  '  Fluid Density           = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity       = ', mu, ' Pa.s'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re)    = ', Re
    write(*,'(A,F12.2,A)')  '  Dynamic Pressure (q)    = ', q_dynamic, ' Pa'
    write(*,*)
    write(*,'(A)') '--- DRAG RESULTS ---'
    write(*,'(A,F12.6)')    '  Drag Coefficient (Cd)   = ', Cd
    write(*,'(A,F12.4,A)')  '  Drag Force (F_D)        = ', F_drag, ' N'
    write(*,'(A,F12.4,A)')  '  Drag Power (P_D)        = ', F_drag * U_inf, ' W'
    write(*,*)
    
    ! Terminal velocity for sphere
    if (shape == 1) then
        write(*,'(A)') '--- TERMINAL VELOCITY (sphere in gravity) ---'
        ! mg = F_drag -> (4/3)pi(D/2)3rho_s g = Cd x 0.5 x rho x V2 x pi(D/2)2
        ! V_t = sqrt(4rho_s g D / (3 Cd rho))
        ! Using rho_s = 2700 kg/m3 (aluminum) as example
        write(*,'(A,F12.4,A)') '  V_terminal (steel)      = ', &
            sqrt(4.0d0 * 7800.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
        write(*,'(A,F12.4,A)') '  V_terminal (aluminum)   = ', &
            sqrt(4.0d0 * 2700.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
        write(*,'(A,F12.4,A)') '  V_terminal (wood)       = ', &
            sqrt(4.0d0 * 600.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
        write(*,*)
    end if
    
    ! Cd vs Re chart data
    write(*,'(A)') '--- Cd vs Re PROFILE ---'
    write(*,'(A)') '  Re            Cd_sphere     Cd_cylinder   Cd_flat'
    write(*,'(A)') '  ---'
    
    do i = 1, 30
        Re_i = 10.0d0**(-1.0d0 + dble(i-1) * 7.0d0 / 29.0d0)
        
        ! Sphere Cd
        if (Re_i < 1.0d0) then
            Cd_i = 24.0d0 / Re_i
        else if (Re_i < 1000.0d0) then
            Cd_i = 24.0d0/Re_i * (1.0d0 + 0.15d0 * Re_i**0.687d0)
        else if (Re_i < 2.0d5) then
            Cd_i = 0.44d0
        else
            Cd_i = 0.1d0
        end if
        
        write(*,'(ES12.4,2X,F12.6,2X,F12.6,2X,F12.6)') Re_i, &
            Cd_i, &
            merge(1.2d0, 1.0d0 + 10.0d0*Re_i**(-2.0d0/3.0d0), Re_i > 1000.0d0), &
            merge(1.17d0, 2.0d0, Re_i > 100.0d0)
    end do
    
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---'
    write(*,'(A)') '  F_D = Cd x (rhoU2/2) x A_ref'
    write(*,'(A)') '  Sphere: Cd varies with Re (Stokes -> Schiller-Naumann)'
    write(*,'(A)') '  Stokes (Re<1): Cd = 24/Re'
    write(*,'(A)') '  Schiller-Naumann: Cd = 24/Re(1 + 0.15 Re^0.687)'
    
end program drag_force_calculator