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