π» Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Lift & Airfoil Calculator
Core Numerical Engine in Fortran 90 β’ 30 total downloads
lift_airfoil.f90
! =========================================================================
! Source File: lift_airfoil.f90
! =========================================================================
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
Solver Description
Analyze NACA 4-digit airfoils using thin airfoil theory with finite-wing corrections. Outputs CL, CD, CM, polar diagram, and airfoil shape.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
π οΈ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 lift_airfoil.f90 -o lift_airfoil
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
lift_airfoil < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! Chord c [m]\nSpan b [m]\nFreestream velocity VΓ’ΛΕΎ [m/s]\nAir density ΓΒ [kg/mΓΒ³]\nViscosity ΓΒΌ [PaΓΒ·s]\nCamber position digit (0Γ’β¬β9)\nnm_init\nnp_init\nnt_init
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0