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
Yes this is possible, and is in fact done all the time. If you want the pole to be at $(x_0,y_0)$ you just make the following substitutions in your equations.
$$ x \rightarrow x-x_0$$ $$ y \rightarrow y-y_0$$
Edit: To answer the question about changing the polar-axis.
One way this can be done is by rotating the coordinates so that the $x-axis$ is at an angle to what it used to be. To rotate the coordinates by and angle $\theta$ apply the transformation,
$$ x' = x \cos(\theta) - y \sin(\theta) $$ $$ y' = x \sin(\theta) + y \sin(\theta) $$
In that case what you would want to substitute into the polar coordinate equations is,
$$ x\rightarrow (x-x_0) \cos(\theta) - (y-y_0) \sin(\theta) $$ $$ y \rightarrow (x-x_0) \sin(\theta) + (y-y_0) \sin(\theta) $$