🔬
Solver Purpose & Physical Scope
Compute Mach angle for supersonic flows. Features interactive shock angle vector visualizations.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 261 times
📄 Source File:
mach_angle.f90
📁 calcul/CFD /
mach_angle.f90
program mach_angle
implicit none
! =================================================================
! Mach Angle / Mach Cone Calculator
! Reads: M, gamma, p1 [Pa], T1 [K], R_gas [J/kg.K], x_dist [m]
! =================================================================
! Constants
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: DEG2RAD = PI / 180.0d0
double precision, parameter :: RAD2DEG = 180.0d0 / PI
! Inputs
double precision :: M, gamma, p1, T1, R_gas, x_dist
! Derived quantities
double precision :: rho1, a1, u1
double precision :: gp1, gm1
! Mach angle
double precision :: mu_rad, mu_deg, full_cone_deg
double precision :: sin_mu, cos_mu, tan_mu
! Cone geometry at distance x
double precision :: cone_radius, cone_diameter, zone_area
! Stagnation conditions
double precision :: T0, p0, rho0
! Disturbance analysis
double precision :: t_sound_1m, dist_source_1m, ratio_dist
! Profile sweep
integer :: i, n_points
double precision :: dM, M_cur, mu_cur, sin_cur, tan_cur, cos_cur
double precision :: p_p0_cur, T_T0_cur, M_sweep_max
integer :: iostat_val
! -----------------------------------------------------------------
! Read inputs
! -----------------------------------------------------------------
read(*,*,iostat=iostat_val) M
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid Mach number input.'
stop
end if
read(*,*,iostat=iostat_val) gamma
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid specific heat ratio input.'
stop
end if
read(*,*,iostat=iostat_val) p1
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid pressure input.'
stop
end if
read(*,*,iostat=iostat_val) T1
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid temperature input.'
stop
end if
read(*,*,iostat=iostat_val) R_gas
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid gas constant input.'
stop
end if
read(*,*,iostat=iostat_val) x_dist
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid distance input.'
stop
end if
! -----------------------------------------------------------------
! Validate inputs
! -----------------------------------------------------------------
if (M < 1.0d0) then
write(*,*) 'ERROR: Mach number must be >= 1.0 for a Mach cone to exist.'
stop
end if
if (gamma <= 1.0d0) then
write(*,*) 'ERROR: Specific heat ratio gamma must be > 1.0.'
stop
end if
if (p1 <= 0.0d0) then
write(*,*) 'ERROR: Pressure must be positive.'
stop
end if
if (T1 <= 0.0d0) then
write(*,*) 'ERROR: Temperature must be positive.'
stop
end if
if (R_gas <= 0.0d0) then
write(*,*) 'ERROR: Gas constant must be positive.'
stop
end if
if (x_dist < 0.0d0) then
write(*,*) 'ERROR: Distance must be >= 0.'
stop
end if
! -----------------------------------------------------------------
! Convenience
! -----------------------------------------------------------------
gp1 = gamma + 1.0d0
gm1 = gamma - 1.0d0
! -----------------------------------------------------------------
! Flow state
! -----------------------------------------------------------------
rho1 = p1 / (R_gas * T1)
a1 = sqrt(gamma * R_gas * T1)
u1 = M * a1
! -----------------------------------------------------------------
! Stagnation conditions
! -----------------------------------------------------------------
T0 = T1 * (1.0d0 + 0.5d0 * gm1 * M**2)
p0 = p1 * (1.0d0 + 0.5d0 * gm1 * M**2) ** (gamma / gm1)
rho0 = rho1 * (1.0d0 + 0.5d0 * gm1 * M**2) ** (1.0d0 / gm1)
! -----------------------------------------------------------------
! Mach angle
! -----------------------------------------------------------------
mu_rad = asin(1.0d0 / M)
mu_deg = mu_rad * RAD2DEG
full_cone_deg = 2.0d0 * mu_deg
sin_mu = sin(mu_rad)
cos_mu = cos(mu_rad)
tan_mu = tan(mu_rad)
! -----------------------------------------------------------------
! Cone geometry at distance x
! -----------------------------------------------------------------
if (x_dist > 0.0d0) then
cone_radius = x_dist * tan_mu
cone_diameter = 2.0d0 * cone_radius
zone_area = PI * cone_radius**2
else
cone_radius = 0.0d0
cone_diameter = 0.0d0
zone_area = 0.0d0
end if
! -----------------------------------------------------------------
! Disturbance analysis
! -----------------------------------------------------------------
t_sound_1m = 1.0d0 / a1 ! time for sound to travel 1 m
dist_source_1m = u1 * t_sound_1m ! distance source moves in that time
ratio_dist = a1 / u1 ! = 1/M = sin(mu)
! -----------------------------------------------------------------
! Print results
! -----------------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' MACH ANGLE / MACH CONE CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- FLOW CONDITIONS -----------------------------------------'
write(*,'(A,F12.4)') ' Mach Number (M) = ', M
write(*,'(A,F12.6)') ' Specific Heat Ratio (g) = ', gamma
write(*,'(A,F12.2,A)') ' Gas Constant (R) = ', R_gas, ' J/kg.K'
write(*,'(A,ES14.6,A)') ' Static Pressure (p) = ', p1, ' Pa'
write(*,'(A,F12.2,A)') ' Static Temperature (T) = ', T1, ' K'
write(*,'(A,F12.4,A)') ' Density (rho) = ', rho1, ' kg/m3'
write(*,'(A,F12.4,A)') ' Speed of Sound (a) = ', a1, ' m/s'
write(*,'(A,F12.4,A)') ' Flow Velocity (u) = ', u1, ' m/s'
write(*,*)
write(*,'(A)') '--- STAGNATION CONDITIONS -----------------------------------'
write(*,'(A,ES14.6,A)') ' Total Pressure p0 = ', p0, ' Pa'
write(*,'(A,F12.2,A)') ' Total Temperature T0 = ', T0, ' K'
write(*,'(A,F12.4,A)') ' Total Density rho0 = ', rho0, ' kg/m3'
write(*,*)
write(*,'(A)') '--- MACH ANGLE RESULTS --------------------------------------'
write(*,'(A,F12.4,A)') ' Mach Angle (mu) = ', mu_deg, ' deg'
write(*,'(A,F12.8,A)') ' Mach Angle (mu) = ', mu_rad, ' rad'
write(*,'(A,F12.4,A)') ' Full Cone Angle (2*mu) = ', full_cone_deg, ' deg'
write(*,'(A,F12.8)') ' sin(mu) = 1/M = ', sin_mu
write(*,'(A,F12.8)') ' cos(mu) = ', cos_mu
write(*,'(A,F12.8)') ' tan(mu) = ', tan_mu
write(*,*)
write(*,'(A)') '--- MACH CONE GEOMETRY --------------------------------------'
write(*,'(A,F12.4,A)') ' Reference Distance (x) = ', x_dist, ' m'
write(*,'(A,F12.6,A)') ' Cone Radius at x = ', cone_radius, ' m'
write(*,'(A,F12.6,A)') ' Cone Diameter at x = ', cone_diameter, ' m'
write(*,'(A,F12.6,A)') ' Zone of Action Area = ', zone_area, ' m2'
write(*,*)
write(*,'(A)') '--- DISTURBANCE PROPAGATION ---------------------------------'
write(*,'(A,ES14.6,A)') ' Time for sound (1 m) = ', t_sound_1m, ' s'
write(*,'(A,F12.4,A)') ' Source travel (in t) = ', dist_source_1m, ' m'
write(*,'(A,F12.6)') ' Ratio a/u = sin(mu) = ', ratio_dist
write(*,*)
! -----------------------------------------------------------------
! Profile sweep: M from 1.0 to M_sweep_max
! -----------------------------------------------------------------
M_sweep_max = max(M * 2.0d0, 10.0d0)
if (M_sweep_max < 2.0d0) M_sweep_max = 10.0d0
write(*,'(A)') '--- MACH ANGLE PROFILE vs MACH NUMBER -----------------------'
write(*,'(A)') ' M mu [deg] sin(mu) tan(mu)' // &
' cos(mu) p/p0 T/T0'
write(*,'(A)') ' ----------------------------------------------------------' // &
'----------------------------'
n_points = 40
dM = (M_sweep_max - 1.0d0) / dble(n_points)
do i = 0, n_points
M_cur = 1.0d0 + dble(i) * dM
if (M_cur < 1.0d0) M_cur = 1.0d0
mu_cur = asin(1.0d0 / M_cur) * RAD2DEG
sin_cur = 1.0d0 / M_cur
cos_cur = sqrt(1.0d0 - sin_cur**2)
if (cos_cur > 1.0d-15) then
tan_cur = sin_cur / cos_cur
else
tan_cur = 0.0d0
end if
T_T0_cur = 1.0d0 / (1.0d0 + 0.5d0 * gm1 * M_cur**2)
p_p0_cur = T_T0_cur ** (gamma / gm1)
write(*,'(F8.4,2X,F11.4,2X,F11.8,2X,F11.8,2X,F11.8,2X,F11.8,2X,F11.8)') &
M_cur, mu_cur, sin_cur, tan_cur, cos_cur, p_p0_cur, T_T0_cur
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' mu = asin(1/M) (Mach angle)'
write(*,'(A)') ' sin(mu) = 1/M = a/u (wave-speed ratio)'
write(*,'(A)') ' tan(mu) = 1/sqrt(M^2-1) (cone radius/distance)'
write(*,'(A)') ' Cone radius at x: r = x * tan(mu)'
write(*,'(A)') ' Full cone angle: 2*mu'
write(*,'(A)') '============================================================'
end program mach_angle
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 mach_angle.f90 -o mach_angle
2. Execution with input.txt redirection:
mach_angle < input.txt
📄 Sample input.txt File Structure
Sample Data:
2.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
Mach Number ($M$)\nSpecific Heat Ratio ($\gamma$)\nStatic Pressure ($p$) [Pa]\nStatic Temperature ($T$) [K]\nGas Preset\nxdist