coordinates - Calculating View Angle?

I am looking for the most accurate way of calculating the view angle (Azimuth, Elevation) from a reference point (my location) A to another point B. I have the Lat, Lon, and Altitude for each point.

Currently I am converting to local ENU coordinates using the Wikipedia formulas for ECEF to ENU and using that to calculate the angles.

Az = atan2(E,N)
El = atan(U / sqrt(E^2 + N^2)

Is there a more accurate way of calculating these Angles?

Maybe a different coordinate system I should be projecting to?

Is this the correct type of Angle, should a rhumb line azimuth or something else be used instead?

I have also looked at using Vincenty's formula for calculating the azimuth but it does not take into account the altitude of each point. A change in Altitude using ENU seems to slightly change the Azimuth which seems correct(?). When comparing Vincenty Azimuth with the ENU Azimuth on the surface (Altitude 0), the answers are extremely close but start to differ if the points are not on the surface.

Answer

For accurate calculations, convert (lat, lon, elevation) directly to earth-centered (x,y,z). (If you don't do this, you need to retain additional information about the local normal ["up"] directions in order to compute angles accurately at nonzero elevations.)

Elevation

Given two points (x,y,z) and (x',y',z') in an earth-centered coordinate system, the vector from the first to the second is (dx,dy,dz) = (x'-x, y'-y, z'-z), whence the cosine of the angle made to the normal at (x,y,z) is the inner product of the unit length versions of those vectors:

Cos(elevation) = (x*dx + y*dy + z*dz) / Sqrt((x^2+y^2+z^2)*(dx^2+dy^2+dz^2))

Obtain its principal inverse cosine. Subtract this from 90 degrees if you want the angle of view relative to a nominal horizon. This is the "elevation."

Azimuth

A similar calculation obtains the local direction of view ("azimuth"). We need a level vector (u,v,w) that points due north. One such vector at the location (x,y,z) is (-zx, -zy, x^2+y^2). (The inner product of these two vectors is zero, proving it is level. Its projection onto the Equatorial plane is proportional to (-x,-y) which points directly inward, making it the projection of a north-pointing vector. These two calculations confirm that this is indeed the desired vector). Therefore

Cos(azimuth) = (-z*x*dx - z*y*dy + (x^2+y^2)*dz) / Sqrt((x^2+y^2)(x^2+y^2+z^2)(dx^2+dy^2+dz^2))

We also need the sine of the azimuth, which is similarly obtained once we know a vector pointing due East (locally). Such a vector is (-y, x, 0), because it clearly is perpendicular to (x,y,z) (the up direction) and the northern direction. Therefore

Sin(azimuth) = (-y*dx + x*dy) / Sqrt((x^2+y^2)(dx^2+dy^2+dz^2))

These values enable us to recover the azimuth as the inverse tangent of the cosine and sine.

---

Example

A pilot in an airplane flying west at 4000 meters, located at (lat, lon) = (39, -75), sees a jet far ahead located at (39, -76) and flying at 12000 meters. What is are the angles of view (relative to the level direction at the pilot's location)?

The XYZ coordinates of the airplanes are (x,y,z) = (1285410, -4797210, 3994830) and (x',y',z') = (1202990, -4824940, 3999870), respectively (in the ITRF00 datum, which uses the GRS80 ellipsoid). The pilot's view vector therefore is (dx,dy,dz) = (-82404.5, -27735.3, 5034.56). Applying the formula gives the cosine of the view angle as 0.0850706. Its inverse cosine is 85.1199 degrees, whence the elevation is 4.88009 degrees: the pilot is looking up by that much.

A north-pointing level vector is (-5.13499, 19.1641, 24.6655) (times 10^12) and an east-pointing level vector is (4.79721, 1.28541, 0) (times 10^6). Therefore, applying the last two formulas, the cosine of the azimuth equals 0.00575112 and its sine equals -0.996358. The ArcTangent function tells us the angle for the direction (0.00575112, -0.996358) is 270.331 degrees: almost due west. (It's not exactly west because the two planes lie on the same circle of latitude, which is curving toward the North pole: see Why is the 'straight line' path across continent so curved? for an extended explanation.)

By the way, this example confirms we got the orientation correct for the azimuth calculation: although it was clear the east-pointing vector was orthogonal to the other two vectors, until now it was not plain that it truly points east and not west.