What is “heading”? To me it means you establish a local coordinate system at point $A_1$, and giving the azimuth (and perhaps elevation) of $A_2$ with respect to this coordinate system. So what directions do you have in this coordinate system?
- Up is radial. The unit up vector is $A_1$, scaled to unit length.
- East is perpendicular to up and perpendicular to the axis of the earth. So compute the cross product of the vector $(0,0,1)$ and the up vector, normalize the result to unit length and make sure you got the correct sign.
- North is perpendicular to up and perpendicular to east, so take their cross product. The norm will already be one (since the other two are length one and perpendicular to one another), so just make sure you got the correct sign.
(Getting all signs correct means taking the cross product in the right order. There are many possible conventions. Instead of discussing this at length, just try out one generic example and if the signs are wrong, swap the order in the cross product.)
So now you have three basis vectors of an orthonormal basis, and you can compute the coordinates of your points with respect to this basis simply by computing the corresponding dot products. If I call these coordinates $u,e,n$ (in direction up, east resp. north), then the azimuth would be
$$\tan\theta=\frac en\qquad\theta=\operatorname{arctan2}(e,n)$$
with zero meaning north and positive angles going east from there. If you point a laser in the given direction at $A_1$ and adjust the elevation correctly, you get it to point at $A_2$ unless something (like the earth) gets in the way.
I interpreted the problem as asking for a function that, for general $R$ and $h$, gives the distance between the two points for any reflection angle (i.e., not only for the case in which the beams are tangent), and allows to calculate the maximal reflection angle and the maximal distance. This function can then be used to calculate the maximal angle and the maximal distance in the specific scenario provided in the OP, with $R=6371$ Km and $h=500$ Km.
Let us consider the Earth circumference as represented by the circle $x^2+y^2=R^2$, with the center in the origin and where $R=6371$. We can place our object in $A(0,R+h)$ on the $y$-axis, representing a point that is $h$ Km up the Earth surface.
Now let us draw two lines passing through $A$, symmetric with respect to the $y$-axis and intersecting the circumference. For each line, let us consider the intersection point that is nearer to the $y$-axis. Let us call the two new points $B$ (in the first quadrant) and $C$ (in the second quadrant). These represent the two points on Earth surface.
Due to the symmetry of the construction, we can continue by analyzing only one of these two points, e.g. $C$. The equation of the line containing $AC$ can be written as $y=sx+R+h$, where $s$ is its positive slope. To determine where this line crosses the circumference, we can set
$$sx+R+h=\sqrt{R^2-x^2}$$
whose solutions are
$$x=\frac{-sR-sh \pm \sqrt{\Delta}}{s^2 + 1}$$
where $\Delta=s^2 R^2 - 2h R - h^2$.
As stated above, we are interested in the less negative solution for $x$, as it is that nearer to the $y$-axis. So we get that the $x$-coordinate of $C$ is
$$X_C=\frac{-sR-sh + \sqrt{\Delta}}{s^2 + 1}$$
and the $y$-coordinate is
$$Y_C=\frac{s(-sR-sh + \sqrt{\Delta})}{s^2 + 1}+R+h$$
As a result, the equation $y=tx$ of the line $OC$ has slope
$$t=\frac{Y_C}{X_C}x=s-\frac{(s^2+1)(R+h)}{sh+sR-\sqrt{\Delta}}$$
Now setting $\angle{BAC}=\alpha$ and $\angle{BOC}=\beta$, we have $s=\cot(\alpha/2)$ and $t=-\cot(\beta/2)$.
Thus, we get
$$ \cot(\beta/2) =\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}}
- \cot(\alpha/2)$$
and then
$$ \beta =2\cot^{-1}\left[\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} - \cot(\alpha/2)\right] $$
So the length of the arc $D$ corresponding to $ \beta$, which is the distance along the spherical surface asked in the OP, is
$$ D =2R\cot^{-1}\left[\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} - \cot(\alpha/2)\right] $$
where $\Delta=\cot^2(\alpha/2)R^2 - 2h R - h^2$.
The last equation can be simplified as
$$ D =2R\cot^{-1}\left[\frac{(R+h)+\sqrt{\Delta}}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} \right] $$
For example, for $\alpha=\pi/2$ and $h= (\sqrt{2}-1)R$, as expected we have $\Delta=0$ (this is the situation where $\alpha$ is a right angle and the light beams are tangent to the surface). In this case, $\beta$ is also a right angle and $D=\pi/2\,R$. Accordingly the formula above gives this result, as shown by WA here.
For any value of $R$ and $h$, the maximal angle $\alpha_{max}$ and the maximal distance $D_{max}$ (i.e., those obtained with the beams tangent to the surface) can be determined by considering the case in which $\Delta=0$. This case occurs when $\cot^2(\alpha/2)R^2 - 2h R - h^2=0$. Solving for $\alpha$ in the range $0 \leq \alpha \leq \pi$ we get
$$\alpha_{max} = 2 \cot^{-1}\left(\frac{\sqrt{h(h+2R)}}{R}\right)$$
Interestingly, when $\Delta=0$, the formula for the distance is considerably simplified, and by few calculations reduces to
$$ D_{max} =2R\cot^{-1}\left[\tan(\alpha/2)\right] $$
As shown here, in the specific scenario described in the OP, substituting $R=6371$ and $h=500$, we get
$$\alpha= 2 \cot^{-1}\left(\frac{10 \sqrt{66210}}{6371}\right) \approx 2.3739 \,\,\text{radians}$$
which corresponds to about $136$ degrees. Here is the plot of the distance $D$ (in Km) as a function of $\alpha$ (in radians) for $R=6371$ and $h=500$, as obtained by WA. The plot confirms the maximal real value of $\alpha$, concordant with the predicted value of $2.3739$. The blue and red lines indicate the real and imaginary part, respectively.
Lastly, from the simplified formula for the maximal distance, taking $R=6371$ and $\alpha=2.3739$, we get
$$D_{max}\approx 4891 \, \text{Km}$$
Best Answer
We measure latitude positive or negative, from the Equator. But I’ll have to talk about the “colatitude”, the distance from the North Pole, so this is $90^\circ-$ ordinary latitude.
Look at your two points, given by colatitude and longitude. Draw geodesics (arcs of great circle) from them to the North Pole, and the geodesic between them. You see that you have an SAS situation (“side-angle-side”), the sides are the lines $\alpha$ and $\beta$ going to the NP, and the angle $\Gamma$ is the difference in the two longitudes. For this, you use the Direct Law of Cosines (there’s a “Polar Law” as well), it goes as follows: $$ \cos\gamma=\cos\alpha\cos\beta+\sin\alpha\sin\beta\cos\Gamma\,, $$ where the lower-case Greek letters are the (angular) measures of sides and the upper-case letters are the angles at the opposite vertices, as is a common convention. That gives your distance $\gamma$ between your two cities, in angular measure of course.
Now, for the angle between two geodesics, you can use the Law of Sines, very easily stated as $$ \frac{\sin\alpha}{\sin A}=\frac{\sin\beta}{\sin B}=\frac{\sin\gamma}{\sin\Gamma}\,. $$ You have the angle $\Gamma$ and now (thanks to your previous computation) the side $\gamma$, so you can very quickly get the angles $A$ and $B$. These can be interpreted as telling you the “initial headings” of a ship sailing from your first city to the second, and vice versa, since they’re angles from due North. When you have the other sides of the triangle, you can add or subtract suitably. You need to draw yourself a picture, of course.
You can also find the angles of your triangle from the three separate distances among the three cities and then use the DLC, solving for the one angle in the formula, but I think LS is easier.