diff --git a/code/hud/hudtarget.cpp b/code/hud/hudtarget.cpp
index 54e14b6b1b6..6d2c7711582 100644
--- a/code/hud/hudtarget.cpp
+++ b/code/hud/hudtarget.cpp
@@ -3135,7 +3135,7 @@ void HudGaugeReticleTriangle::renderTriangle(vec3d *hostile_pos, int aspect_flag
 		unsize(&hostile_vertex.screen.xyw.x, &hostile_vertex.screen.xyw.y);
 	}
 
-	float ang = atan2(-(hostile_vertex.screen.xyw.y - tablePosY), hostile_vertex.screen.xyw.x - tablePosX);
+	float ang = atan2_safe(-(hostile_vertex.screen.xyw.y - tablePosY), hostile_vertex.screen.xyw.x - tablePosX);
 	float sin_ang=sinf(ang);
 	float cos_ang=cosf(ang);
 
diff --git a/code/math/vecmat.cpp b/code/math/vecmat.cpp
index dc878e33f9a..80aec695b7a 100644
--- a/code/math/vecmat.cpp
+++ b/code/math/vecmat.cpp
@@ -63,6 +63,20 @@ bool vm_matrix_equal(const matrix4 &self, const matrix4 &other)
 		vm_vec_equal(self.vec.pos, other.vec.pos);
 }
 
+// -----------------------------------------------------------
+// atan2_safe()
+//
+// Wrapper around atan2() that handles the special case of x == 0 and y == 0.
+//
+float atan2_safe(float y, float x)
+{
+	// special case
+	if ( x == 0.0f && y == 0.0f )
+		return 0.0f;
+
+	return atan2(y, x);
+}
+
 // ---------------------------------------------------------------------
 // vm_vec_component()
 //
@@ -1109,7 +1123,7 @@ angles *vm_extract_angles_matrix(angles *a, const matrix *m)
 {
 	float sinh,cosh,cosp;
 
-	a->h = atan2(m->vec.fvec.xyz.x,m->vec.fvec.xyz.z);
+	a->h = atan2_safe(m->vec.fvec.xyz.x,m->vec.fvec.xyz.z);
 
 	sinh = sinf(a->h); cosh = cosf(a->h);
 
@@ -1124,7 +1138,7 @@ angles *vm_extract_angles_matrix(angles *a, const matrix *m)
 
 	fvec_xz_distance = fl_sqrt( ( (m->vec.fvec.xyz.x)*(m->vec.fvec.xyz.x) ) + ( (m->vec.fvec.xyz.z)*(m->vec.fvec.xyz.z) ) );
 
-	a->p = atan2(-m->vec.fvec.xyz.y, fvec_xz_distance);
+	a->p = atan2_safe(-m->vec.fvec.xyz.y, fvec_xz_distance);
 
 	if (cosp == 0.0f)	//the cosine of pitch is zero.  we're pitched straight up. say no bank
 
@@ -1136,7 +1150,7 @@ angles *vm_extract_angles_matrix(angles *a, const matrix *m)
 		sinb = m->vec.rvec.xyz.y/cosp;
 		cosb = m->vec.uvec.xyz.y/cosp;
 
-		a->b = atan2(sinb,cosb);
+		a->b = atan2_safe(sinb,cosb);
 	}
 
 
@@ -1182,7 +1196,7 @@ static angles *vm_extract_angles_vector_normalized(angles *a, const vec3d *v)
 
 	a->p = asinf_safe(-v->xyz.y);
 
-	a->h = atan2(v->xyz.z,v->xyz.x);
+	a->h = atan2_safe(v->xyz.z,v->xyz.x);
 
 	return a;
 }
@@ -1853,10 +1867,10 @@ float vm_closest_angle_to_matrix(const matrix* mat, const vec3d* rot_axis, float
 		//If we support IEEE float handling, we don't need this, the div by 0 will be handled correctly with the INF. If not, do this:
 		const float yz_recip = (!std::numeric_limits<float>::is_iec559 && y * z < 0.001f) ? FLT_MAX : 1.0f / (y * z);
 
-		solutions = { 2.0f * atan2f(-sr_neg * (y * y + sr) * yz_recip, -2.0f * sr_neg),
-					  2.0f * atan2f(sr_neg * (y * y + sr) * yz_recip, 2.0f * sr_neg),
-					  2.0f * atan2f(-sr_pos * (y * y - sr) * yz_recip, -2.0f * sr_pos),
-					  2.0f * atan2f(sr_pos * (y * y - sr) * yz_recip, 2.0f * sr_pos) };
+		solutions = { 2.0f * atan2_safe(-sr_neg * (y * y + sr) * yz_recip, -2.0f * sr_neg),
+					  2.0f * atan2_safe(sr_neg * (y * y + sr) * yz_recip, 2.0f * sr_neg),
+					  2.0f * atan2_safe(-sr_pos * (y * y - sr) * yz_recip, -2.0f * sr_pos),
+					  2.0f * atan2_safe(sr_pos * (y * y - sr) * yz_recip, 2.0f * sr_pos) };
 	}
 	float value = -2.0f;
 	float correct = 0;
diff --git a/code/math/vecmat.h b/code/math/vecmat.h
index f3fe9c2958b..cc5cc7cb6dc 100755
--- a/code/math/vecmat.h
+++ b/code/math/vecmat.h
@@ -509,6 +509,9 @@ void vm_angular_move_forward_vec(const vec3d *goal_fvec, const matrix *orient, c
 // Find the bounding sphere for a set of points (center and radius are output parameters)
 void vm_find_bounding_sphere(const vec3d *pnts, int num_pnts, vec3d *center, float *radius);
 
+// Wrapper around atan2() that returns 0 for the degenerate (0, 0) input.
+float atan2_safe(float y, float x);
+
 // Translates from world coordinates to body coordinates
 vec3d* vm_rotate_vec_to_body(vec3d *body_vec, const vec3d *world_vec, const matrix *orient);
 
diff --git a/code/render/3ddraw.cpp b/code/render/3ddraw.cpp
index 08364688a31..eb1b4ab373e 100644
--- a/code/render/3ddraw.cpp
+++ b/code/render/3ddraw.cpp
@@ -893,7 +893,7 @@ void g3_render_laser_2d(material *mat_params, vec3d *headp, float head_width, ve
 		w = len_2d;
 
 	} else {
-		a = atan2(taily - heady, tailx - headx);
+		a = atan2_safe(taily - heady, tailx - headx);
 
 		w = len_2d;
 
diff --git a/code/render/3dmath.cpp b/code/render/3dmath.cpp
index 609d7a55980..9ed8a6b12a2 100644
--- a/code/render/3dmath.cpp
+++ b/code/render/3dmath.cpp
@@ -94,10 +94,10 @@ ubyte g3_transfer_vertex(vertex *dest, const vec3d *src)
 static void g3_compensate_asymmetric_fov(float& x, float& y, float z) {
 	if (std::holds_alternative<asymmetric_fov>(Proj_fov)) {
 		const auto& afov = std::get<asymmetric_fov>(Proj_fov);
-		float angle = atan2(z, x) + (afov.left + afov.right);
+		float angle = atan2_safe(z, x) + (afov.left + afov.right);
 		x = angle == PI_2 ? 0.0f : z / tanf(angle);
 
-		angle = atan2(z, y) + (afov.up + afov.down);
+		angle = atan2_safe(z, y) + (afov.up + afov.down);
 		y = angle == PI_2 ? 0.0f : z / tanf(angle);
 	}
 }