Ignore the points $A,B,C,D$ in the following sketch.
The translation corresponds to the change of variables
$$x=\overline{x}+X,\qquad y=\overline{y}+Y,\tag{1}$$
and the rotation corresponds to
$$X=x^{\prime }\cos \theta -y^{\prime }\sin \theta,\qquad Y=x^{\prime }\sin \theta +y^{\prime }\cos \theta .\tag{2}$$
If you combine $(1)$ and $(2)$, you get
$$x=\overline{x}+x^{\prime }\cos \theta -y^{\prime }\sin \theta,\qquad y=\overline{y}+x^{\prime }\sin \theta +y^{\prime }\cos \theta .\tag{3}$$
You can invert this system of equations to get $x',y'$ in terms of $x,y$
and $\overline{x},\overline{y}$:
$$x^{\prime }=x\cos \theta +y\sin \theta -\overline{x}\cos \theta -\overline{y}\sin \theta \tag{4a}$$
$$y^{\prime }=-x\sin \theta +y\cos \theta +\overline{x}\sin \theta -\overline{y}\cos
\theta .\tag{4b}$$
The rotation and the translation transformation can be expressed by matrices. The rotation transformation in two dimensions is described in this Wikipedia section.
I haven't checked the formula, but you have ${850000 \over 6371000}$ instead of
either ${8500000 \over 6371000}$ or ${8500 \over 6371}$.
The point for the -180 bearing in the data above is about 750 km south of the center location.
(I believe the correct coordinates are about $(-96.1953, -75.7800)$.)
It is not clear to me how you are computing your new coordinates. The Haversine formula is used to compute the distance given two sets of coordinates and does
not have a trivial inverse. Perhaps you can provide a link?
Here is one way to do it, I am sure there must be better ways:
I am performing the computations in the x-y-z plane rather than using
latitude/longitude. The mapping between the two is straightforward.
Given a point $p \in \mathbb{R}^3$ on the surface of the earth excluding the poles, a bearing $\beta$ (measured clockwise from north) and a distance $L$,
compute a new point $p'$ a distance $L$ away from $p$ at bearing $\beta$ (from $p$).
Note that the poles are excluded so that the bearing makes sense, at either pole a bearing to/from North makes no sense.
I am assuming a perfectly spherical Earth. To reduce notational clutter,
define
$u: \mathbb{R}^3 \setminus \{0\} \to \mathbb{R}^3$ by $u(x) = {1 \over \|x\|} x$ (that is, normalise $x$).
Let $R= \|p\|$, be the radius (of the idealised Earth), and let $N=(0,0,R)$ be the porth pole.
Note that a distance $L$ corresponds to an angle (at the Earth's center) of
$\alpha = {L \over R}$ (radians).
Compute a 'local north direction' $\nu = u(N-{1 \over R^2}\langle p, N \rangle p)$. This is a direction contained in the intersection of the tangent plane at $p$ with the subspace containing $p$ and $N$.
Compute a 'local east direction' $e = u(\nu \times p)$. This is a
direction perpendicular to both $\nu$ and $e$ with the appropriate
orientation.
(Note that $\nu, p, e$ are mutually orthogonal.)
Compute a 'local direction' $d = (\sin \beta) e + (\cos \beta) \nu$.
This is a unit vector in the direction we want to go in on the tangent
plane at $p$.
Then we have the new point $p' = (\cos \alpha) p + (\sin \alpha) R d$.
Best Answer
Assuming that your bearing is measured east of north (right of the $y$-axis), you can use $$ \text{bearing}=2\;\tan^{-1}\left(\frac{x}{y+\sqrt{x^2+y^2}}\right)\tag{1} $$ $\text{bearing}<0^\circ$ means west of north.
$\text{bearing}>90^\circ$ means south of east.
$\text{bearing}<-90^\circ$ means south of west.
Formula $(1)$ only fails when $y\le0$ and $x=0$ (division by $0$); in that case, $\text{bearing}=180^\circ$, due south.