Accuracy of GPS bearing, given two locations and their accuracies For an Android app I need to determine the accuracy of the bearing returned by the GPS. Android supplies the following data, among others:


*

*The GPS position

*The accuracy of the GPS position, expressed in meters, which corresponds to $\sigma$, assuming a Gaussian distribution

*The bearing


Note that no accuracy information is supplied for the bearing. Since the only thing the GPS receiver calculates directly is the position, while bearing is inferred from two subsequent position, my approach is to take two subsequent positions and their accuracy to determine an accuracy for the bearing.
Let


*

*$P_0$ and $P_1$ be two subsequent positions,

*$\epsilon_0$ and $\epsilon_1$ their respective accuracies in meters,

*$\psi_1$ the bearing reported at $P_1$,

*$d$ the distance between $P_0$ and $P_1$ in meters,

*$\psi_\epsilon$ the unknown accuracy of the bearing, expressed in the same unit as the bearing (degrees or radians at your option)


We can visualize this as a line segment from $P_0$ to $P_1$ having length $d$, and circles around each point, their radii being $\epsilon_0$ and $\epsilon_1$, respectively.
As long as $\epsilon_0 + \epsilon_1 <= d$, $\psi_\epsilon$ reaches its maximum when the actual locations are located just on these circles, so that the line segment connecting them lies on the common tangent of the two circles which intersects $d$.
Together with the two line segments connecting each measured location with the (presumed) actual location, this becomes a figure of two similar rectangular triangles.
Then:
$$\epsilon_0 + \epsilon_1 = d * \sin (\psi_\epsilon)$$
and thus:
$$\psi_\epsilon = \arcsin (\frac{\epsilon_0 + \epsilon_1}{d})$$
The limit is reached when $\epsilon_0 + \epsilon_1 = d$, in which case $\psi_\epsilon$ is 90°.
As soon as $\epsilon_0 + \epsilon_1 > d$, the above formula is no longer defined. Visualizing the situation, the accuracy circles overlap, which means that the actual locations can be any pair of points that is entirely within both circles and thus any bearing is possible, hence $\psi_\epsilon = 180°$.
The issues here is the sudden "jump" from 90° to 180°; also I probably haven't paid proper attention to probability here.
Is there a better formula for calculating bearing accuracy, given two positions and their accuracy?
EDIT:
It's been pointed out that bearing accuracy is influenced by several factors. The one that springs to my mind is that the GPS assumes the receiver has travelled along a straight line from $P_0$ and $P_1$, which will result in errors while turning. Those need to be considered separately – for now, I just need to know how the accuracy of the bearing depends on positional accuracy, without challenging any of the assumptions made by the GPS.
 A: I don't like this problem because the accuracy in bearing given by the GPS depends on may factors and is calculated using also doppler effect on the received signals (so is usually very accurate, as long as the speed is large enough).
Moreover, to calculate the variances on positions and bearings it is better to use tools as Kalman filters (which GPSs already use).
Anyway, to give an answer following your directions I think you can use the definition of differential.
Your bearing can be calculated as: $$ b=\tan^{-1}\left ( \frac{\Delta y}{\Delta x} \right ) $$
Its differential is: $$db=\frac{\partial b}{\partial x}dx+\frac{\partial b}{\partial y}dy=-\frac{\Delta y}{d^2}dx+\frac{\Delta x}{d^2}dy$$ 
and represents a linear approximation of the variation of the bearing function function as its inputs change.
If we assume $\epsilon_1+\epsilon_2$ as the variation of the inputs: $$ \psi_e = \left ( \frac{\Delta x-\Delta y}{d^2} \right ) (\epsilon_1+\epsilon_2)$$ Where $\psi_e$ is in radians. This function can assume negative values so you can just take its modulo to have always a positive error. 
Anyway, if this worries you, we can use the four quadrant version of the inverse tangent: $$ b=2\tan^{-1}\left (\frac{d-\Delta x}{\Delta y} \right )$$
obtaining:$$ \psi_e = \left ( \frac{\Delta x+\Delta y - d}{d(d-\Delta x)} \right ) (\epsilon_1+\epsilon_2)$$
As $\Delta x+\Delta y \gt d$ and $\Delta x \lt d$ this equation is positively definite.
