💻 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.
Oblique Shock Wave relations
Core Numerical Engine in Fortran 90 • 59 total downloads
oblique_shock.f90
! =========================================================================
! Source File: oblique_shock.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc — Oblique Shock Wave Relations Solver
! Reference: NACA Report 1135, Anderson Modern Compressible Flow Ch. 4
! ============================================================================
program oblique_shock
implicit none
! Input variables
integer :: solve_for ! 1=Solve for beta (Weak), 2=Solve for beta (Strong), 3=Solve for theta, 4=Solve for M1
double precision :: M1 ! Upstream Mach number
double precision :: theta ! Deflection angle [deg]
double precision :: beta ! Shock angle [deg]
double precision :: gamma ! Specific heat ratio (default 1.4)
! Constants
double precision, parameter :: pi = 3.141592653589793d0
double precision, parameter :: rad = pi / 180.0d0
! Intermediate and output variables
double precision :: theta_rad, beta_rad, mu, mu_deg
double precision :: beta_star_rad, beta_star_deg, theta_max_rad, theta_max_deg
double precision :: Mn1, Mn2, M2
double precision :: P2_P1, T2_T1, rho2_rho1, P02_P01
logical :: detached
! Solver variables
double precision :: low, high, mid, g_mid, val
integer :: iter
! Read input from stdin
read(*,*) solve_for
read(*,*) M1
read(*,*) theta
read(*,*) beta
read(*,*) gamma
! Check inputs
if (gamma <= 1.0d0) gamma = 1.4d0
detached = .false.
! Initialize defaults for output ratios
P2_P1 = 1.0d0
T2_T1 = 1.0d0
rho2_rho1 = 1.0d0
P02_P01 = 1.0d0
M2 = 0.0d0
Mn1 = 0.0d0
Mn2 = 0.0d0
mu = 0.0d0
mu_deg = 0.0d0
beta_star_rad = 0.0d0
beta_star_deg = 0.0d0
theta_max_deg = 0.0d0
! ── CASE 1, 2, 3: Upstream Mach checks ───────────────────
if (solve_for /= 4) then
if (M1 <= 1.0d0) then
write(*,*) "ERROR: Upstream Mach number M1 must be supersonic (> 1.0) for oblique shocks."
stop
end if
! Calculate Mach Angle
mu = asin(1.0d0 / M1)
mu_deg = mu / rad
! Calculate beta* and theta_max analytically
val = (gamma + 1.0d0) * (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2 + (gamma + 1.0d0)/16.0d0 * M1**4)
val = (gamma + 1.0d0)/4.0d0 * M1**2 - 1.0d0 + sqrt(val)
val = val / (gamma * M1**2)
if (val >= 0.0d0 .and. val <= 1.0d0) then
beta_star_rad = asin(sqrt(val))
beta_star_deg = beta_star_rad / rad
theta_max_rad = atan(f_beta(beta_star_rad, M1, gamma))
theta_max_deg = theta_max_rad / rad
end if
end if
! ── SOLVE ACCORDING TO MODE ──────────────────────────────
select case (solve_for)
case (1, 2) ! Solve for beta (Weak or Strong)
theta_rad = theta * rad
if (theta > theta_max_deg) then
detached = .true.
beta_rad = beta_star_rad ! Default to max deflection limit
beta = beta_star_deg
else
if (solve_for == 1) then
! Weak solution on [mu, beta_star]
low = mu
high = beta_star_rad
do iter = 1, 150
mid = 0.5d0 * (low + high)
g_mid = atan(f_beta(mid, M1, gamma)) - theta_rad
if (g_mid < 0.0d0) then
low = mid
else
high = mid
end if
if (abs(high - low) < 1.0d-14) exit
end do
beta_rad = mid
beta = beta_rad / rad
else
! Strong solution on [beta_star, 90 deg]
low = beta_star_rad
high = 90.0d0 * rad
do iter = 1, 150
mid = 0.5d0 * (low + high)
g_mid = atan(f_beta(mid, M1, gamma)) - theta_rad
if (g_mid > 0.0d0) then
low = mid
else
high = mid
end if
if (abs(high - low) < 1.0d-14) exit
end do
beta_rad = mid
beta = beta_rad / rad
end if
end if
case (3) ! Solve for theta given M1 and beta
beta_rad = beta * rad
if (beta_rad <= mu) then
write(*,*) "ERROR: Shock angle beta must be greater than the Mach angle mu = ", mu_deg
stop
end if
if (beta >= 90.0d0) then
theta = 0.0d0
theta_rad = 0.0d0
else
val = f_beta(beta_rad, M1, gamma)
if (val < 0.0d0) then
write(*,*) "ERROR: Given shock angle results in non-physical negative deflection."
stop
end if
theta_rad = atan(val)
theta = theta_rad / rad
end if
case (4) ! Solve for M1 given theta and beta
theta_rad = theta * rad
beta_rad = beta * rad
if (beta <= theta) then
write(*,*) "ERROR: Shock angle beta must be greater than deflection angle theta."
stop
end if
! Solve for M1^2 using rearranged equation
val = sin(2.0d0 * beta_rad) - tan(theta_rad) * (gamma + cos(2.0d0 * beta_rad))
if (val <= 0.0d0) then
write(*,*) "ERROR: No physical Mach number solution exists for the given angles (denominator is non-positive)."
stop
end if
mid = 2.0d0 * (1.0d0 / tan(beta_rad) + tan(theta_rad)) / val
if (mid <= 1.0d0) then
write(*,*) "ERROR: Calculated Mach number M1 must be supersonic (> 1.0)."
stop
end if
M1 = sqrt(mid)
mu = asin(1.0d0 / M1)
mu_deg = mu / rad
! Calculate beta* and theta_max for solved M1
val = (gamma + 1.0d0) * (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2 + (gamma + 1.0d0)/16.0d0 * M1**4)
val = (gamma + 1.0d0)/4.0d0 * M1**2 - 1.0d0 + sqrt(val)
val = val / (gamma * M1**2)
if (val >= 0.0d0 .and. val <= 1.0d0) then
beta_star_rad = asin(sqrt(val))
beta_star_deg = beta_star_rad / rad
theta_max_rad = atan(f_beta(beta_star_rad, M1, gamma))
theta_max_deg = theta_max_rad / rad
end if
case default
write(*,*) "ERROR: Invalid solve target mode."
stop
end select
! ── COMPUTE DOWNSTREAM PROPERTIES (IF ATTACHED) ──────────
if (.not. detached) then
beta_rad = beta * rad
theta_rad = theta * rad
Mn1 = M1 * sin(beta_rad)
if (Mn1 > 1.0d0) then
! Normal Shock Relations applied to Mn1
Mn2 = sqrt((2.0d0 + (gamma - 1.0d0) * Mn1**2) / (2.0d0 * gamma * Mn1**2 - (gamma - 1.0d0)))
M2 = Mn2 / sin(beta_rad - theta_rad)
P2_P1 = (2.0d0 * gamma * Mn1**2 - (gamma - 1.0d0)) / (gamma + 1.0d0)
rho2_rho1 = ((gamma + 1.0d0) * Mn1**2) / ((gamma - 1.0d0) * Mn1**2 + 2.0d0)
T2_T1 = P2_P1 / rho2_rho1
P02_P01 = (rho2_rho1)**(gamma / (gamma - 1.0d0)) * (1.0d0 / P2_P1)**(1.0d0 / (gamma - 1.0d0))
else
! Sonic/Subsonic normal component (should not occur if inputs are physical)
M2 = M1
P2_P1 = 1.0d0
rho2_rho1 = 1.0d0
T2_T1 = 1.0d0
P02_P01 = 1.0d0
end if
end if
! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
write(*, '(A, A)') "Detachment Status = ", merge("Detached", "Attached", detached)
write(*, '(A, F14.6)') "Shock Angle (beta) = ", beta
write(*, '(A, F14.6)') "Deflection Angle (theta) = ", theta
write(*, '(A, F14.6)') "Upstream Mach (M1) = ", M1
if (detached) then
write(*, '(A, F14.6)') "Downstream Mach (M2) = ", 0.0d0
write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", 0.0d0
write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", 0.0d0
write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", 0.0d0
write(*, '(A, F14.6)') "Total Pressure Ratio (P02/P01) = ", 0.0d0
else
write(*, '(A, F14.6)') "Downstream Mach (M2) = ", M2
write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
write(*, '(A, F14.6)') "Total Pressure Ratio (P02/P01) = ", P02_P01
end if
write(*, '(A, F14.6)') "Max Deflection Angle (theta_max) = ", theta_max_deg
write(*, '(A, F14.6)') "Shock Angle Max Deflection (beta_max) = ", beta_star_deg
write(*, '(A, F14.6)') "Mach Angle (mu) = ", mu_deg
write(*, '(A, F14.6)') "Upstream Normal Mach (Mn1) = ", Mn1
write(*, '(A, F14.6)') "Downstream Normal Mach (Mn2) = ", Mn2
write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma
contains
! Helper function to evaluate the theta-beta-Mach relation: tan(theta)
double precision function f_beta(b, M, g)
double precision, intent(in) :: b ! Shock angle [rad]
double precision, intent(in) :: M ! Mach number
double precision, intent(in) :: g ! Specific heat ratio
f_beta = 2.0d0 * (1.0d0 / tan(b)) * (M**2 * (sin(b))**2 - 1.0d0) / &
(M**2 * (g + cos(2.0d0 * b)) + 2.0d0)
end function f_beta
end program oblique_shock
Solver Description
Solve oblique shock wave relations (theta-beta-Mach relation). Supports attached/detached checks, weak/strong solutions, and interactive curves.
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 oblique_shock.f90 -o oblique_shock_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
oblique_shock_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Solve target (1=beta weak, 2=beta strong, 3=theta, 4=M1)
1
! Upstream Mach number M1
2.0
! Deflection angle theta [deg]
10.0
! Shock wave angle beta [deg]
40.0
! Specific heat ratio gamma
1.4
1
! Upstream Mach number M1
2.0
! Deflection angle theta [deg]
10.0
! Shock wave angle beta [deg]
40.0
! Specific heat ratio gamma
1.4