💻 Fortran Source Code Library

! =========================================================================
! Source File: optimal_ld.f90
! =========================================================================

!==============================================================================
! ThermoFluidCalc — Calculator #25 : Optimal L/D Ratio
!==============================================================================
! Physics : For a parabolic drag polar  CD = CD0 + CL^2 / (pi*e*AR) :
!
!   Max Range   (max L/D)   :  CDi = CD0
!     CL = sqrt(pi*e*AR*CD0)   ,   (L/D)_max = 1/(2*sqrt(CD0/(pi*e*AR)))
!
!   Max Endurance (min CL^(3/2)/CD)  :  CDi = 3*CD0
!     CL = sqrt(3*pi*e*AR*CD0)
!
!   Max Climb Rate (min excess power) :  CDi = CD0/3
!     CL = sqrt(pi*e*AR*CD0/3)
!
!   V_opt = sqrt(2*W / (rho*S*CL_opt))
!
! Reference : Gupta, Appendix 7.2
!
! Build:
!   gfortran -O2 -o optimal_ld optimal_ld.f90
!
! Input (stdin, one line):
!   CD0  e  AR  rho  S  W  CL_max  npts_sensitivity
!==============================================================================
program optimal_ld
  implicit none

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

  integer  :: npts, i
  real(dp) :: CD0, e, AR, rho, S, W, CL_max
  real(dp) :: k   ! k = 1/(pi*e*AR)
  ! Max range
  real(dp) :: CL_r, CD_r, LD_r, V_r
  ! Max endurance
  real(dp) :: CL_e, CD_e, LD_e, V_e
  ! Max climb
  real(dp) :: CL_c, CD_c, LD_c, V_c
  ! Min speed
  real(dp) :: V_min, CD_min, LD_min
  ! Sensitivity
  real(dp) :: cd0_lo, cd0_hi, dcd0, cd0_v, cl_v, ld_v

  read(*,*) CD0, e, AR, rho, S, W, CL_max, npts

  ! Validate
  if (CD0 <= 0.0_dp) then; write(*,'(A)') 'ERROR=CD0 must be positive.'; stop; end if
  if (e <= 0.0_dp .or. e > 1.0_dp) then; write(*,'(A)') 'ERROR=Oswald e must be (0,1].'; stop; end if
  if (AR <= 0.0_dp) then; write(*,'(A)') 'ERROR=AR must be positive.'; stop; end if
  if (rho <= 0.0_dp) then; write(*,'(A)') 'ERROR=Density must be positive.'; stop; end if
  if (S <= 0.0_dp) then; write(*,'(A)') 'ERROR=Wing area must be positive.'; stop; end if
  if (W <= 0.0_dp) then; write(*,'(A)') 'ERROR=Weight must be positive.'; stop; end if
  if (CL_max <= 0.0_dp) then; write(*,'(A)') 'ERROR=CL_max must be positive.'; stop; end if
  if (npts < 5) npts = 5
  if (npts > MAX_PTS) npts = MAX_PTS

  k = 1.0_dp / (PI * e * AR)

  ! ── 1. Max Range (CDi = CD0) ──────────────────────────────────────────
  CL_r = sqrt(CD0 / k)       ! = sqrt(pi*e*AR*CD0)
  CD_r = 2.0_dp * CD0
  LD_r = CL_r / CD_r
  V_r  = sqrt(2.0_dp * W / (rho * S * CL_r))

  ! ── 2. Max Endurance (CDi = 3*CD0) ────────────────────────────────────
  CL_e = sqrt(3.0_dp * CD0 / k)
  CD_e = 4.0_dp * CD0
  LD_e = CL_e / CD_e
  V_e  = sqrt(2.0_dp * W / (rho * S * CL_e))

  ! ── 3. Max Climb Rate (CDi = CD0/3) ──────────────────────────────────
  CL_c = sqrt(CD0 / (3.0_dp * k))
  CD_c = CD0 + CD0 / 3.0_dp   ! = 4/3 CD0
  LD_c = CL_c / CD_c
  V_c  = sqrt(2.0_dp * W / (rho * S * CL_c))

  ! ── 4. Minimum speed ──────────────────────────────────────────────────
  V_min  = sqrt(2.0_dp * W / (rho * S * CL_max))
  CD_min = CD0 + CL_max**2 * k
  if (CD_min > 0.0_dp) then
    LD_min = CL_max / CD_min
  else
    LD_min = 0.0_dp
  end if

  ! ── Output ────────────────────────────────────────────────────────────
  write(*,'(A,ES15.8)') 'CD0=', CD0
  write(*,'(A,F8.4)')   'OSWALD_E=', e
  write(*,'(A,F10.4)')  'AR=', AR
  write(*,'(A,ES15.8)') 'RHO=', rho
  write(*,'(A,ES15.8)') 'S=', S
  write(*,'(A,ES15.8)') 'W=', W
  write(*,'(A,F8.4)')   'CL_MAX=', CL_max
  write(*,'(A,ES15.8)') 'K=', k

  ! Range
  write(*,'(A,F10.6)')  'CL_RANGE=', CL_r
  write(*,'(A,F10.6)')  'CD_RANGE=', CD_r
  write(*,'(A,F10.4)')  'LD_RANGE=', LD_r
  write(*,'(A,F10.3)')  'V_RANGE=', V_r

  ! Endurance
  write(*,'(A,F10.6)')  'CL_ENDUR=', CL_e
  write(*,'(A,F10.6)')  'CD_ENDUR=', CD_e
  write(*,'(A,F10.4)')  'LD_ENDUR=', LD_e
  write(*,'(A,F10.3)')  'V_ENDUR=', V_e

  ! Climb
  write(*,'(A,F10.6)')  'CL_CLIMB=', CL_c
  write(*,'(A,F10.6)')  'CD_CLIMB=', CD_c
  write(*,'(A,F10.4)')  'LD_CLIMB=', LD_c
  write(*,'(A,F10.3)')  'V_CLIMB=', V_c

  ! Min speed
  write(*,'(A,F10.3)')  'V_MIN=', V_min
  write(*,'(A,F10.6)')  'CD_VMIN=', CD_min
  write(*,'(A,F10.4)')  'LD_VMIN=', LD_min

  ! ── Sensitivity sweep ─────────────────────────────────────────────────
  cd0_lo = 0.5_dp * CD0
  cd0_hi = 2.0_dp * CD0
  dcd0   = (cd0_hi - cd0_lo) / real(npts - 1, dp)

  write(*,'(A)') 'SENSITIVITY_START'
  do i = 0, npts - 1
    cd0_v = cd0_lo + real(i, dp) * dcd0
    cl_v  = sqrt(cd0_v / k)
    ld_v  = cl_v / (2.0_dp * cd0_v)
    write(*,'(ES15.8,A,F12.4)') cd0_v, ',', ld_v
  end do
  write(*,'(A)') 'SENSITIVITY_END'

end program optimal_ld