program boundary_layer
    implicit none
    
    double precision :: U_inf, rho, nu, L, Re_cr
    double precision :: Re_L, x_tr, Re_x
    double precision :: delta, delta_star, theta, Cf, tau_w
    double precision :: q_inf, CD_total, F_drag, A_corr
    double precision :: pi, x, dx
    integer :: i, n_stations
    character(len=20) :: regime
    
    ! Turbulent virtual origin
    double precision :: x_virt
    double precision :: delta_lam_tr, delta_turb
    
    pi = 3.14159265358979d0
    n_stations = 40
    
    ! Read inputs
    read(*,*) U_inf      ! Free-stream velocity [m/s]
    read(*,*) rho        ! Fluid density [kg/m3]
    read(*,*) nu         ! Kinematic viscosity [m2/s]
    read(*,*) L          ! Plate length [m]
    read(*,*) Re_cr      ! Critical Reynolds number for transition
    
    ! Global calculations
    Re_L = U_inf * L / nu
    x_tr = Re_cr * nu / U_inf
    q_inf = 0.5d0 * rho * U_inf**2
    
    ! ============================================
    ! HEADER
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   BOUNDARY LAYER ANALYSIS — FLAT PLATE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ---'
    write(*,'(A,F12.4,A)')  '  Free-stream velocity U∞ = ', U_inf, ' m/s'
    write(*,'(A,F12.4,A)')  '  Fluid density ρ         = ', rho, ' kg/m³'
    write(*,'(A,ES12.4,A)') '  Kinematic viscosity ν   = ', nu, ' m²/s'
    write(*,'(A,F12.4,A)')  '  Plate length L          = ', L, ' m'
    write(*,'(A,ES12.4)')   '  Critical Re_x           = ', Re_cr
    write(*,*)
    write(*,'(A)') '--- GLOBAL RESULTS ---'
    write(*,'(A,ES12.4)')   '  Reynolds at plate end Re_L = ', Re_L
    
    if (x_tr < L) then
        write(*,'(A,F12.6,A)') '  Transition location x_tr   = ', x_tr, ' m'
        write(*,'(A,F8.2,A)')  '  Transition at              = ', (x_tr/L)*100.0d0, ' % of plate'
        write(*,'(A)')         '  Flow type                  = Mixed (laminar + turbulent)'
    else
        write(*,'(A,F12.6,A)') '  Transition location x_tr   = ', x_tr, ' m (beyond plate)'
        write(*,'(A)')         '  Flow type                  = Fully laminar'
    end if
    write(*,*)
    
    ! ============================================
    ! TOTAL DRAG COEFFICIENT
    ! ============================================
    write(*,'(A)') '--- TOTAL DRAG COEFFICIENT ---'
    
    if (Re_L < Re_cr) then
        ! Fully laminar
        CD_total = 1.328d0 / sqrt(Re_L)
        write(*,'(A,ES12.6)')   '  C_D (laminar)           = ', CD_total
        write(*,'(A)')          '  Formula: C_D = 1.328 / √Re_L'
    else
        ! Mixed boundary layer
        ! A depends on Re_cr: A = Re_cr × (0.074/Re_cr^0.2 − 1.328/√Re_cr)
        A_corr = Re_cr * (0.074d0 / Re_cr**0.2d0 - 1.328d0 / sqrt(Re_cr))
        CD_total = 0.074d0 / Re_L**0.2d0 - A_corr / Re_L
        write(*,'(A,ES12.6)')   '  C_D (mixed BL)          = ', CD_total
        write(*,'(A,F12.4)')    '  Correction A            = ', A_corr
        write(*,'(A)')          '  Formula: C_D = 0.074/Re_L^0.2 − A/Re_L'
    end if
    
    F_drag = CD_total * q_inf * L * 1.0d0  ! per unit width
    write(*,'(A,F12.6,A)')  '  Drag force (per m width)= ', F_drag, ' N/m'
    write(*,'(A,F12.4,A)')  '  Dynamic pressure q∞     = ', q_inf, ' Pa'
    write(*,*)
    
    ! ============================================
    ! BOUNDARY LAYER AT PLATE END
    ! ============================================
    write(*,'(A)') '--- BOUNDARY LAYER AT x = L ---'
    Re_x = U_inf * L / nu
    
    if (Re_L < Re_cr) then
        delta = 5.0d0 * L / sqrt(Re_x)
        delta_star = 1.7208d0 * L / sqrt(Re_x)
        theta = 0.664d0 * L / sqrt(Re_x)
        Cf = 0.664d0 / sqrt(Re_x)
        regime = 'Laminar'
    else
        delta = 0.37d0 * L / Re_x**0.2d0
        delta_star = delta / 8.0d0
        theta = 7.0d0 * delta / 72.0d0
        Cf = 0.0592d0 / Re_x**0.2d0
        regime = 'Turbulent'
    end if
    
    tau_w = Cf * q_inf
    
    write(*,'(A,A)')        '  Regime                  = ', trim(regime)
    write(*,'(A,F12.6,A)')  '  δ  (BL thickness)       = ', delta*1000, ' mm'
    write(*,'(A,F12.6,A)')  '  δ* (displacement)       = ', delta_star*1000, ' mm'
    write(*,'(A,F12.6,A)')  '  θ  (momentum)           = ', theta*1000, ' mm'
    write(*,'(A,ES12.6)')   '  C_f (skin friction)     = ', Cf
    write(*,'(A,F12.4,A)')  '  τ_w (wall shear)        = ', tau_w, ' Pa'
    write(*,'(A,F12.6)')    '  Shape factor H = δ*/θ   = ', delta_star / theta
    write(*,*)
    
    ! ============================================
    ! PROFILE TABLE
    ! ============================================
    write(*,'(A)') '--- BOUNDARY LAYER PROFILE ALONG PLATE ---'
    write(*,'(A)') '  x [m]      Re_x         δ [mm]     δ* [mm]    θ [mm]     Cf          τ_w [Pa]   Regime'
    write(*,'(A)') '  -------    ----------   --------   --------   --------   ----------  --------   ------'
    
    do i = 1, n_stations
        x = L * dble(i) / dble(n_stations)
        Re_x = U_inf * x / nu
        
        if (Re_x < Re_cr) then
            ! Laminar (Blasius)
            delta = 5.0d0 * x / sqrt(Re_x)
            delta_star = 1.7208d0 * x / sqrt(Re_x)
            theta = 0.664d0 * x / sqrt(Re_x)
            Cf = 0.664d0 / sqrt(Re_x)
            regime = 'Laminar'
        else
            ! Turbulent (1/7 power law)
            delta = 0.37d0 * x / Re_x**0.2d0
            delta_star = delta / 8.0d0
            theta = 7.0d0 * delta / 72.0d0
            Cf = 0.0592d0 / Re_x**0.2d0
            regime = 'Turbulent'
        end if
        
        tau_w = Cf * q_inf
        
        write(*,'(F8.4,2X,ES12.4,2X,F8.4,3X,F8.4,3X,F8.4,3X,ES10.4,2X,F8.4,3X,A)') &
            x, Re_x, delta*1000, delta_star*1000, theta*1000, Cf, tau_w, trim(regime)
    end do
    
    write(*,*)
    
    ! ============================================
    ! CHART DATA (machine-readable)
    ! ============================================
    write(*,'(A)') '--- CHART_DATA ---'
    do i = 1, n_stations
        x = L * dble(i) / dble(n_stations)
        Re_x = U_inf * x / nu
        
        if (Re_x < Re_cr) then
            delta = 5.0d0 * x / sqrt(Re_x)
            delta_star = 1.7208d0 * x / sqrt(Re_x)
            theta = 0.664d0 * x / sqrt(Re_x)
            Cf = 0.664d0 / sqrt(Re_x)
            write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,ES12.6,A,I1)') &
                x, ',', delta*1000, ',', delta_star*1000, ',', theta*1000, ',', Cf, ',', 0
        else
            delta = 0.37d0 * x / Re_x**0.2d0
            delta_star = delta / 8.0d0
            theta = 7.0d0 * delta / 72.0d0
            Cf = 0.0592d0 / Re_x**0.2d0
            write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,ES12.6,A,I1)') &
                x, ',', delta*1000, ',', delta_star*1000, ',', theta*1000, ',', Cf, ',', 1
        end if
    end do
    write(*,'(A)') '--- END_CHART_DATA ---'
    
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  LAMINAR (Blasius):'
    write(*,'(A)') '    δ = 5.0 x / √Re_x'
    write(*,'(A)') '    δ* = 1.7208 x / √Re_x'
    write(*,'(A)') '    θ = 0.664 x / √Re_x'
    write(*,'(A)') '    C_f = 0.664 / √Re_x'
    write(*,'(A)') ''
    write(*,'(A)') '  TURBULENT (1/7 power law):'
    write(*,'(A)') '    δ = 0.37 x / Re_x^0.2'
    write(*,'(A)') '    δ* = δ / 8'
    write(*,'(A)') '    θ = 7δ / 72'
    write(*,'(A)') '    C_f = 0.0592 / Re_x^0.2'
    write(*,'(A)') ''
    write(*,'(A)') '  MIXED DRAG:'
    write(*,'(A)') '    C_D = 0.074/Re_L^0.2 − A/Re_L'
    write(*,'(A)') '    where A = Re_cr(0.074/Re_cr^0.2 − 1.328/√Re_cr)'
    
end program boundary_layer