Is my perspective math correct?
Identifying your axes in both figures and adding the camera position to your first figure would help you understand what's going on.

You could also have a single variables for all your points, generating a 2D matrix with the rows as each point and columns as the components $x$, $y$ and $z$. That way, you could handle the projection using a simple matrix multiplication instead of handling each row separately.

In the problem statement, would be interesting to know the direction of your rotation, and more importantly, the original camera direction and its up vector. I guess you camera is rotated by 160° counter-clockwise, the original camera direction is $[0, 0, 1]$ and its up vector is $[0, 1, 0]$. If any of those assumptions is wrong, the rest of the answer will be wrong.

There is no sign of mm to m conversion in your code. Either your focal should be $0.016$ and $S_x = S_y = 0.0001$ or your point coordinates should be multiplied by $0.00001$.

Let's reason about where should the points end up in your image. For example, your camera center, without rotation, is pointing at the $[-1,1,x]$ line. Since the points are all on the $z=0.5$ plane, we can do the following analysis: if we focus on the $x$-axis, we can see that $tan(160°) \cdot (5 - 0.5) = 1.64...$, so the camera center will end up slightly to the left of the points (as the camera is at $x=-1$), so the center will end up at $\approx 0.64$, meaning that the points will appear on the right section of the image. Furthermore, the camera has the same $y$ coordinates as two point, and since the $y$ coordinates are not changed by the rotation, they should still end up on the same coordinates after the transform, meaning on the central row of the image.

A good way to check your answer is by using an existing 3D modeler like Blender:

Be cautious with Blender' coordinates system, for example the default camera vector is [0, 0, -1].
Here's the render:

Focal was set to another value to make the sphere more visible. So we see that the bottom two points are on the middle row of the image and the points are slightly to the right of the image.

I implemented your homework in Python:

import numpy as np

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D


# Parameters
f_mm = 0.016
f_px = f_mm / 0.00001
t_cam = np.array([[-1., 1., 5.]]).T
t_cam_homogeneous = np.vstack((t_cam, np.array([[0]])))
theta = 160. * np.pi / 180.
ox = 320
oy = 240
# Rotation and points are in homogeneous coordinates
rot_cam = np.array([[np.cos(theta), 0, np.sin(theta)],
                    [0, 1, 0],
                    [-np.sin(theta), 0, np.cos(theta)]])
points = np.array([[1, 1, 0.5, 1],
                   [1, 1.5, 0.5, 1],
                   [1.5, 1.5, 0.5, 1],
                   [1.5, 1, 0.5, 1]]).T

# Compute projection matrix using intrinsics and extrinsics
intrinsics = np.array([[f_px, 0, ox],
                       [0, f_px, oy],
                       [0, 0, 1]])
extrinsics = np.hstack((rot_cam, rot_cam.dot(-t_cam)))

rot_cam2 = np.identity(4); rot_cam2[:3,:3] = rot_cam
camera_coordinates = rot_cam2.dot(points - t_cam_homogeneous)
camera_coordinates = camera_coordinates[:3,:] / camera_coordinates[3,:]

# Perform the projection
projected_points = intrinsics.dot(camera_coordinates)
projected_points = projected_points[:2,:] / projected_points[2,:]
projected_points[0,:] = -projected_points[0,:] # Inverted x-axis because camera is pointing toward [0, 0, 1]

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(points[0,:], points[1,:], points[2,:], label="Points")
ax.scatter(t_cam[0], t_cam[1], t_cam[2], c="red", label="Camera")
ax.set_xlabel("X axis"); ax.set_ylabel("Y axis"); ax.set_zlabel("Z axis")
plt.title("World coordinates")
plt.legend()
plt.savefig('world_coordinates.png', dpi=300, bbox_inches="tight")

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(camera_coordinates[0,:], camera_coordinates[1,:], camera_coordinates[2,:], label="Points")
ax.scatter(0, 0, 0, c="red", label="Camera")
ax.set_xlabel("X axis"); ax.set_ylabel("Y axis"); ax.set_zlabel("Z axis")
plt.title("Camera coordinates")
plt.legend()
plt.savefig('camera_coordinates.png', dpi=300, bbox_inches="tight")

plt.figure()
plt.scatter(projected_points[0,:], projected_points[1,:])
plt.xlabel("X axis"); plt.ylabel("Y axis")
plt.title("Image coordinates")
plt.savefig('image_coordinates.png', dpi=300, bbox_inches="tight")

plt.show()


This yield me those figures:




Respectively: World coordinates, Camera coordinates, Camera coordinates rotated to fit slightly the camera orientation (note that here the camera vector goes toward the figure viewpoint, it does not "enter" the figure) and image coordinates.

So we see that the vertical coordinates for the bottom points are correctly on the middle row (240) and the points are on the right side of the image (horizontal value > 320).

I believe one bug you had is that you found negative X values, so you negated the focals (-f/Sxy) in the intrinsics matrix to compensate. The issue here is that we've assumed that the camera initially was pointing toward $[0, 0, 1]$ (otherwise the 160° rotation would not point toward the points). If you look at it that way, the $x$-axis increases when going toward the left, the inverse of this axis should be taken.

Both our results seems similar to me, just that you supposed a $[0, -1, 0]$ up vector for the camera (in fact both axes were mirrored because you negated both focals), and did the computation in mm instead of meters.