!==============================================================================
! ThermoFluidCalc — Calculator #24 : Lift CL & Induced Drag CDi
!==============================================================================
! Physics :
!   CL  = CL_alpha * (alpha - alpha_0)       (linear lift model)
!   CDi = CL^2 / (pi * e * AR)               (induced drag)
!   CD  = CD0 + CDi                           (total drag, parabolic polar)
!   AR  = b^2 / S                             (aspect ratio)
!   L   = CL * q * S ,   D = CD * q * S
!   L/D = CL / CD
!
! Modes:
!   1 = Single point (one angle of attack)
!   2 = Alpha sweep  (range of angles)
!
! Reference : Gupta, §7.8, Eq. 7.87
!
! Build:
!   gfortran -O2 -o lift_drag lift_drag.f90
!==============================================================================
program lift_drag
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: PI = 3.141592653589793238_dp
  integer, parameter :: MAX_PTS = 5000

  integer  :: mode, npts, i
  real(dp) :: rho, V, S, b, e, CD0, CL_alpha, alpha0_deg, alpha_deg
  real(dp) :: alpha_min, alpha_max, dalpha
  real(dp) :: AR, q, CL, CDi, CD, L_force, D_force, LD
  real(dp) :: alpha_rad, alpha0_rad
  real(dp) :: best_LD, best_alpha, best_CL

  read(*,*) mode

  select case (mode)

  !==========================================================================
  ! MODE 1 : Single point
  !==========================================================================
  case (1)
    backspace(5)
    read(*,*) mode, rho, V, S, b, e, CD0, CL_alpha, alpha_deg, alpha0_deg

    call validate_inputs(rho, V, S, b, e)

    AR        = b**2 / S
    q         = 0.5_dp * rho * V**2
    alpha_rad = alpha_deg * PI / 180.0_dp
    alpha0_rad= alpha0_deg * PI / 180.0_dp
    CL        = CL_alpha * (alpha_rad - alpha0_rad)
    CDi       = CL**2 / (PI * e * AR)
    CD        = CD0 + CDi
    L_force   = CL * q * S
    D_force   = CD * q * S
    if (CD > 0.0_dp) then
      LD = CL / CD
    else
      LD = 0.0_dp
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Single Point'
    write(*,'(A,F12.6)')    'AR=', AR
    write(*,'(A,ES15.8)')   'Q=', q
    write(*,'(A,F12.6)')    'ALPHA_DEG=', alpha_deg
    write(*,'(A,F12.6)')    'ALPHA0_DEG=', alpha0_deg
    write(*,'(A,F12.6)')    'CL=', CL
    write(*,'(A,F12.6)')    'CDI=', CDi
    write(*,'(A,F12.6)')    'CD0=', CD0
    write(*,'(A,F12.6)')    'CD=', CD
    write(*,'(A,ES15.8)')   'LIFT=', L_force
    write(*,'(A,ES15.8)')   'DRAG=', D_force
    write(*,'(A,F12.4)')    'LD=', LD

  !==========================================================================
  ! MODE 2 : Alpha sweep
  !==========================================================================
  case (2)
    backspace(5)
    read(*,*) mode, rho, V, S, b, e, CD0, CL_alpha, alpha0_deg, &
              alpha_min, alpha_max, npts

    call validate_inputs(rho, V, S, b, e)
    if (npts < 2) npts = 2
    if (npts > MAX_PTS) npts = MAX_PTS

    AR         = b**2 / S
    q          = 0.5_dp * rho * V**2
    alpha0_rad = alpha0_deg * PI / 180.0_dp
    dalpha     = (alpha_max - alpha_min) / real(npts - 1, dp)

    best_LD    = -1.0e30_dp
    best_alpha = alpha_min
    best_CL    = 0.0_dp

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Alpha Sweep'
    write(*,'(A,F12.6)')    'AR=', AR
    write(*,'(A,ES15.8)')   'Q=', q
    write(*,'(A,F12.6)')    'CD0=', CD0
    write(*,'(A,I5)')       'NPTS=', npts

    write(*,'(A)') 'DATA_START'
    do i = 0, npts - 1
      alpha_deg = alpha_min + real(i, dp) * dalpha
      alpha_rad = alpha_deg * PI / 180.0_dp
      CL  = CL_alpha * (alpha_rad - alpha0_rad)
      CDi = CL**2 / (PI * e * AR)
      CD  = CD0 + CDi
      if (CD > 0.0_dp) then
        LD = CL / CD
      else
        LD = 0.0_dp
      end if

      if (LD > best_LD) then
        best_LD    = LD
        best_alpha = alpha_deg
        best_CL    = CL
      end if

      write(*,'(F8.3,4(A,F12.6))') alpha_deg, ',', CL, ',', CDi, ',', CD, ',', LD
    end do
    write(*,'(A)') 'DATA_END'

    write(*,'(A,F12.4)')   'MAX_LD=', best_LD
    write(*,'(A,F8.3)')    'ALPHA_MAXLD=', best_alpha
    write(*,'(A,F12.6)')   'CL_MAXLD=', best_CL

  case default
    write(*,'(A)') 'ERROR=Invalid mode (must be 1 or 2).'
    stop
  end select

contains

  subroutine validate_inputs(rho, V, S, b, e)
    real(dp), intent(in) :: rho, V, S, b, e
    if (rho <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Density must be positive.'; stop
    end if
    if (V <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Velocity must be positive.'; stop
    end if
    if (S <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Wing area must be positive.'; stop
    end if
    if (b <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Wingspan must be positive.'; stop
    end if
    if (e <= 0.0_dp .or. e > 1.0_dp) then
      write(*,'(A)') 'ERROR=Oswald efficiency must be in (0, 1].'; stop
    end if
  end subroutine

end program lift_drag