Airfoil Vortex Panel & Induced Drag

program airfoil_panel_method_vortex
    implicit none
    integer :: iostat_val
    double precision :: chord_mm, AR_aspect, oswald_e, alpha_deg, m_camber_pct, p_pos_tenth, t_thick_pct
    double precision :: alpha_rad, alpha0_deg, alpha0_rad, a0_slope, a3D_slope
    double precision :: CL_val, CD0_val, CDi_val, CD_tot, L_D_ratio, CM_c4
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) chord_mm       ! Airfoil Chord [mm] (e.g. 1000.0)
    read(*,*,iostat=iostat_val) AR_aspect      ! Wing Aspect Ratio AR (e.g. 8.0)
    read(*,*,iostat=iostat_val) oswald_e       ! Oswald Efficiency Factor e (e.g. 0.85)
    read(*,*,iostat=iostat_val) alpha_deg      ! Angle of Attack [deg] (e.g. 5.0)
    read(*,*,iostat=iostat_val) m_camber_pct   ! NACA Max Camber [% chord] (e.g. 2.0)
    read(*,*,iostat=iostat_val) p_pos_tenth    ! NACA Camber Position (e.g. 4.0)
    read(*,*,iostat=iostat_val) t_thick_pct    ! NACA Max Thickness [% chord] (e.g. 12.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for airfoil panel vortex calculation.'
        stop
    end if

    if (chord_mm <= 0.0d0 .or. AR_aspect <= 0.5d0 .or. oswald_e <= 0.0d0) then
        write(*,*) 'ERROR: Chord, aspect ratio, and Oswald factor must be positive.'
        stop
    end if

    alpha_rad = alpha_deg * (PI / 180.0d0)

    ! Zero-lift angle of attack for NACA 4-digit camber line [deg]
    alpha0_deg = -1.15d0 * m_camber_pct
    alpha0_rad = alpha0_deg * (PI / 180.0d0)

    ! 2D theoretical lift slope a0 = 2*pi [1/rad]
    a0_slope = 2.0d0 * PI * (1.0d0 + 0.77d0 * (t_thick_pct / 100.0d0))

    ! 3D Finite Wing Lift Slope (Prandtl Lifting Line Theory)
    a3D_slope = a0_slope / (1.0d0 + (a0_slope / (PI * AR_aspect)))

    ! Total Lift Coefficient CL
    CL_val = a3D_slope * (alpha_rad - alpha0_rad)

    ! Parasitic Drag CD0 (Hoerner skin friction + form drag)
    CD0_val = 0.0065d0 * (1.0d0 + 2.0d0 * (t_thick_pct / 100.0d0) + 60.0d0 * ((t_thick_pct / 100.0d0)**4))

    ! Induced Drag CDi
    CDi_val = (CL_val**2) / (PI * oswald_e * AR_aspect)
    CD_tot = CD0_val + CDi_val

    if (CD_tot > 0.0001d0) then
        L_D_ratio = CL_val / CD_tot
    else
        L_D_ratio = 0.0d0
    end if

    ! Quarter-chord pitching moment CM_c/4
    CM_c4 = -(PI / 2.0d0) * (m_camber_pct / 100.0d0) * (1.0d0 - p_pos_tenth / 10.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — AIRFOIL VORTEX PANEL & LIFT/DRAG ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2)') 'Chord / Aspect Ratio AR   = ', chord_mm, ' mm / ', AR_aspect
    write(*,'(A,F10.1,A,F10.1,A)') 'NACA Profile / Camber     = ', m_camber_pct, ' % / ', p_pos_tenth*10.0d0, ' % chord'
    write(*,'(A,F10.2,A,F10.2,A)') 'Angle of Attack / Zero-L  = ', alpha_deg, ' deg / ', alpha0_deg, ' deg'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4)')   'LIFT COEFFICIENT (CL)     = ', CL_val
    write(*,'(A,F10.4)')   'TOTAL DRAG COEFF (CD)     = ', CD_tot
    write(*,'(A,F10.4,A,F10.4,A)') 'Parasitic / Induced Drag  = ', CD0_val, ' (CD0) / ', CDi_val, ' (CDi)'
    write(*,'(A,F10.2)')   'LIFT-TO-DRAG RATIO (L/D)  = ', L_D_ratio
    write(*,'(A,F10.4)')   'Moment Coeff (CM_c/4)     = ', CM_c4
    write(*,'(A,F10.4,A)') '3D Lift Slope (dCL/dalpha)= ', a3D_slope * (PI / 180.0d0), ' per degree'
    write(*,'(A)') '============================================================'

end program airfoil_panel_method_vortex