program lift_airfoil
    implicit none
    integer :: i, n_aoa, iostat_val, naca_m, naca_p, naca_t
    double precision :: chord, span, Vinf, rho, mu, alpha_deg, alpha_rad
    double precision :: CL, CD, CM, L, D, M_moment, Re_c, AR, q_inf
    double precision :: CL0, CLa, alpha_L0, alpha_stall, CD0, CD_ind, e_oswald
    double precision :: a_i, a_rad, CL_i, CD_i, CM_i, L_i, D_i
    double precision :: m_cam, p_cam, t_thick
    double precision, parameter :: PI = 3.141592653589793d0
    character(len=20) :: naca_str

    read(*,*,iostat=iostat_val) chord
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid chord input.'
        stop
    end if
    read(*,*,iostat=iostat_val) span
    read(*,*,iostat=iostat_val) Vinf
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) alpha_deg
    read(*,*,iostat=iostat_val) naca_m
    read(*,*,iostat=iostat_val) naca_p
    read(*,*,iostat=iostat_val) naca_t
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all airfoil inputs.'
        stop
    end if
    if (chord <= 0.0d0 .or. span <= 0.0d0 .or. Vinf <= 0.0d0) then
        write(*,*) 'ERROR: Chord, span, 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

    ! NACA 4-digit interpretation
    m_cam = dble(naca_m) / 100.0d0      ! max camber fraction
    p_cam = dble(naca_p) / 10.0d0       ! position of max camber
    t_thick = dble(naca_t) / 100.0d0    ! max thickness fraction
    write(naca_str,'(I1,I1,I2.2)') naca_m, naca_p, naca_t

    alpha_rad = alpha_deg * PI / 180.0d0
    Re_c = rho * Vinf * chord / mu
    AR = span / chord
    q_inf = 0.5d0 * rho * Vinf**2

    ! Thin airfoil theory: CL_alpha = 2*pi per radian
    CLa = 2.0d0 * PI
    ! Zero-lift angle from camber
    if (p_cam > 0.0d0 .and. m_cam > 0.0d0) then
        alpha_L0 = -2.0d0 * m_cam * (1.0d0 - 2.0d0*p_cam)
    else
        alpha_L0 = 0.0d0
    end if
    CL0 = CLa * (-alpha_L0)

    ! Lift coefficient (2D)
    CL = CLa * (alpha_rad - alpha_L0)

    ! 3D finite-wing correction (lifting line)
    if (AR > 0.5d0) then
        CL = CL / (1.0d0 + CLa/(PI*AR))
    end if

    ! Stall estimate (simplified)
    alpha_stall = 12.0d0 + 2.0d0*m_cam*100.0d0
    if (alpha_stall > 20.0d0) alpha_stall = 20.0d0
    if (abs(alpha_deg) > alpha_stall) then
        ! Post-stall: CL drops, use Viterna flat-plate-like
        CL = CL * (alpha_stall / max(abs(alpha_deg), 0.1d0))
    end if

    ! Drag coefficient
    ! Parasitic drag from flat plate estimate
    if (Re_c > 500.0d0) then
        CD0 = 0.074d0 / Re_c**0.2d0
    else
        CD0 = 1.328d0 / sqrt(max(Re_c, 1.0d0))
    end if
    ! Form factor for thickness
    CD0 = CD0 * (1.0d0 + 2.0d0*t_thick + 60.0d0*t_thick**4)
    ! Induced drag
    e_oswald = 0.85d0
    if (AR > 0.5d0) then
        CD_ind = CL**2 / (PI * e_oswald * AR)
    else
        CD_ind = 0.0d0
    end if
    CD = CD0 + CD_ind

    ! Post-stall drag increase
    if (abs(alpha_deg) > alpha_stall) then
        CD = CD + 0.5d0 * sin(2.0d0*alpha_rad)**2
    end if

    ! Moment coefficient about quarter-chord (thin airfoil)
    if (p_cam > 0.0d0 .and. m_cam > 0.0d0) then
        CM = -PI/2.0d0 * m_cam * (1.0d0 - 2.0d0*p_cam)
    else
        CM = 0.0d0
    end if

    ! Forces
    L = CL * q_inf * chord * span
    D = CD * q_inf * chord * span
    M_moment = CM * q_inf * chord**2 * span

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   LIFT & AIRFOIL CALCULATOR ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- AIRFOIL DATA --------------------------------------------'
    write(*,'(A,A)')        '  NACA Profile             = ', trim(naca_str)
    write(*,'(A,F10.4)')    '  Max Camber (m)           = ', m_cam
    write(*,'(A,F10.4)')    '  Camber Position (p)      = ', p_cam
    write(*,'(A,F10.4)')    '  Max Thickness (t/c)      = ', t_thick
    write(*,'(A,ES12.4,A)') '  Chord                    = ', chord, ' m'
    write(*,'(A,ES12.4,A)') '  Span                     = ', span, ' m'
    write(*,'(A,ES12.4)')   '  Aspect Ratio             = ', AR
    write(*,*)
    write(*,'(A)') '--- OPERATING CONDITIONS ------------------------------------'
    write(*,'(A,ES12.4,A)') '  Freestream Velocity      = ', Vinf, ' m/s'
    write(*,'(A,ES12.4,A)') '  Density                  = ', rho, ' kg/m3'
    write(*,'(A,ES12.4)')   '  Reynolds Number          = ', Re_c
    write(*,'(A,F10.3,A)')  '  Angle of Attack          = ', alpha_deg, ' deg'
    write(*,'(A,F10.3,A)')  '  Zero-Lift Angle          = ', alpha_L0*180.0d0/PI, ' deg'
    write(*,'(A,F10.3,A)')  '  Stall Angle Estimate     = ', alpha_stall, ' deg'
    write(*,*)
    write(*,'(A)') '--- AERODYNAMIC COEFFICIENTS --------------------------------'
    write(*,'(A,ES12.4)')   '  Lift Coefficient CL      = ', CL
    write(*,'(A,ES12.4)')   '  Drag Coefficient CD      = ', CD
    write(*,'(A,ES12.4)')   '  Parasitic Drag CD0       = ', CD0
    write(*,'(A,ES12.4)')   '  Induced Drag CDi         = ', CD_ind
    write(*,'(A,ES12.4)')   '  Moment Coefficient CM    = ', CM
    write(*,'(A,ES12.4)')   '  Oswald Efficiency        = ', e_oswald
    write(*,'(A,ES12.4)')   '  L/D Ratio                = ', CL/max(CD,1.0d-30)
    write(*,*)
    write(*,'(A)') '--- FORCES AND MOMENT ---------------------------------------'
    write(*,'(A,ES12.4,A)') '  Lift Force               = ', L, ' N'
    write(*,'(A,ES12.4,A)') '  Drag Force               = ', D, ' N'
    write(*,'(A,ES12.4,A)') '  Pitching Moment          = ', M_moment, ' N.m'
    write(*,*)

    ! Polar sweep: AoA from -8 to +22 degrees
    n_aoa = 61
    write(*,'(A)') '--- POLAR DATA ----------------------------------------------'
    write(*,'(A)') '  AoA[deg]     CL            CD            CM            L/D'
    write(*,'(A)') '  ----------------------------------------------------------------'
    do i = 1, n_aoa
        a_i = -8.0d0 + 30.0d0*dble(i-1)/dble(n_aoa-1)
        a_rad = a_i * PI / 180.0d0

        CL_i = CLa * (a_rad - alpha_L0)
        if (AR > 0.5d0) CL_i = CL_i / (1.0d0 + CLa/(PI*AR))
        if (abs(a_i) > alpha_stall) &
            CL_i = CL_i * (alpha_stall / max(abs(a_i), 0.1d0))

        if (AR > 0.5d0) then
            CD_i = CD0 + CL_i**2 / (PI*e_oswald*AR)
        else
            CD_i = CD0
        end if
        if (abs(a_i) > alpha_stall) &
            CD_i = CD_i + 0.5d0*sin(2.0d0*a_rad)**2

        CM_i = CM

        write(*,'(F10.3,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            a_i, CL_i, CD_i, CM_i, CL_i/max(CD_i,1.0d-30)
    end do
    write(*,*)

    ! Airfoil shape coordinates
    write(*,'(A)') '--- AIRFOIL SHAPE -------------------------------------------'
    write(*,'(A)') '  x/c           y_upper/c     y_lower/c'
    write(*,'(A)') '  -------------------------------------------'
    do i = 0, 50
        call naca4_coords(dble(i)/50.0d0, m_cam, p_cam, t_thick)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Thin airfoil: CL = 2*pi*(alpha - alpha_L0).'
    write(*,'(A)') '  Finite wing: CL_3D = CL_2D / (1 + 2*pi/(pi*AR)).'
    write(*,'(A)') '  Induced drag: CDi = CL^2 / (pi e AR).'
    write(*,'(A)') '  NACA 4-digit thickness and camber distributions.'

contains

    subroutine naca4_coords(xc, m, p, tc)
        implicit none
        double precision, intent(in) :: xc, m, p, tc
        double precision :: yt, yc_cam, dyc, theta, xu, yu, xl, yl
        ! Thickness distribution
        yt = 5.0d0*tc*(0.2969d0*sqrt(max(xc,0.0d0)) &
             - 0.1260d0*xc - 0.3516d0*xc**2 &
             + 0.2843d0*xc**3 - 0.1015d0*xc**4)
        ! Camber line
        if (p > 0.0d0 .and. m > 0.0d0) then
            if (xc < p) then
                yc_cam = m/(p**2)*(2.0d0*p*xc - xc**2)
                dyc = 2.0d0*m/(p**2)*(p - xc)
            else
                yc_cam = m/((1.0d0-p)**2)*((1.0d0-2.0d0*p) &
                         + 2.0d0*p*xc - xc**2)
                dyc = 2.0d0*m/((1.0d0-p)**2)*(p - xc)
            end if
        else
            yc_cam = 0.0d0
            dyc = 0.0d0
        end if
        theta = atan(dyc)
        xu = xc - yt*sin(theta)
        yu = yc_cam + yt*cos(theta)
        xl = xc + yt*sin(theta)
        yl = yc_cam - yt*cos(theta)
        write(*,'(F10.6,2X,F12.6,2X,F12.6)') xc, yu, yl
    end subroutine naca4_coords

end program lift_airfoil